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FOREWORD 


This is the final r eport descr ib ing a stu dy conducted for the National Aero- 
nautic s and^pace^Atliiunlstration , George C. Marshall Space Flight Center 
under Contract Number NAS 8-28607. This study was a "Parachute 
Dynamics and Stability Analysis" as applied to the Solid Rocket Booster 
recovery system of the Space Shuttle. This report covers the period from 
1 February 1973 through 1 February 1974. The Contract Technical Monitor 
is Mr. Gaines L. Watts 

The authors wish to express their gratitude to Mr. M. Bazakos for his assis- 
tance with the computer simulation programs that were developed and to 
Dr. R. E. Rose, Program Manager, for his guidance and supervision. 
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PARACHUTE DYNAMICS AND STABILITY ANALYSIS 


By: S. K. Ibrahim and R. A. Engdahl 


SUMMARY 


The nonlinear differential equations of motion for a general parachute-riser- 
payload system are developed. The resulting math model is then applied for 
analyzing the descent dynamics and stability characteristics ot both the 
drogue stabilization phase and the main descent phase of the Space Shuttle 
Solid Rocket Booster (SRB) recovery system. 

The formulation of the problem is characterized by a minimum number of 
simplifying assumptions and full application of state-of-the-art parachute 
technology. The parachute suspension lines and the parachute risers can be 
modeled as elastic elements, and the whole system may be subjected to 
specified wind and gust profiles in order to assess their effects on the sta- 
bility of the recovery system. 

A numerical linearization technique is provided as an optional subroutine. 

It permits the linearization of the system's equations of motion at selected 
points of the descent trajectory and the calculation of the Eigenvalues 
describing the principal motions. Root locus plots may be obtained to study 
the variation in stability characteristics as a function of time. Computer 
simulations with the nonlinear system of equations were run for a wide range 
of initial conditions both with and without the elastic suspension system 
effects and the wind and gust models. For selected runs, the linearization 
procedure was applied at predetermined points, the Eigenvalues were cal- 
culated, and the stability characteristics were examined. It was determined 
that, for the range of anticipated initial conditions, the projected drogue 
configuration quickly stabilizes the SRB motions, the SRB/ Main descent con- 
figuration is stable, and the motions of the system, with the specified wind 
and gust profiles, remain within acceptable limits at water impact. 


INTRODUCTION 


This is the final report of a one -year program of analytical and computa- 
tional work. The program's primary objective was to formulate a realistic 
mathematical model for the descent dynamics of a parachute /vehicle system 
and to use that model as the basis for a computer simulation, stability 
analysis, and parametric optimization of the Space Shuttle Solid Rocket 
Booster (SRB) recovery system. 
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The recoverable weight of the SRB is at least three times that of any pre- 
viously recovered payload. Full scale testing may not be feasible and large 
scale drop tests are very costly; hence, the need for realistic simulation 
models to permit detailed studies of optimum system parameters and sta- 
bility characteristics and to minimize the number of drop tests. 

The math model described in this report is more general than previously 
published models. Among other things, it permits 6 degrees -of -freedom 
motion for both the parachtue and the vehicle, it includes elastic represen- 
tation for the risers and suspension line, the effect of deterministic winds 
and gusts on the system's performance and a more general representation 
of apparent mass effects. A separate computer program, using the elastic 
element approach, permits the calculation of more realistic canopy profile 
shapes. 


LIST OF SYMBOLS 


ALCM 

B -t 

lk 


B S1" B S3 

B S4" B S6 

C. 

1 

C Ni 

CM 

r 

^M t 


CP 


F li 


3i 


TV 


g 

I. 


Length from confluence point to plane of skirt 

Direction cosines matrix element, body j, row i, column 
k (i, j, k = 1, 2, 3) 

Direction cosine scalar products. Parachute 
Direction cosine scalar products, SRB 
Velocity vector of body i, i = 1, 2, 3 
Normal force coefficient, body i, i ~ 1, 3 
Center of mass 

Moment Coefficient, body i, i = 1, 3 
Tangent force coefficient, body i, i = 1, 3 
Center of pressure 

Aerodynamic forces on body 1 in direction i(i=X, Y, Z), lb 
Aerodynamic forces on body 3 in direction i(i=X, Y, Z), lb 
Riser force, lb 

2 

Gravitational acceleration, ft /sec 

o 

Principal moments of inertia matrix, body i.slug ft 
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I 


A1 
K 


LS 


K 

L 

L, 


R 


L V L R 


J 3T 


"4 


J CM 


m i 


m T 


m lA 


m. 


m. 


M 3i 

N 


P i 

Qi 

R. 


R 


2 

Principal apparent moments of inertia matrix, slug ft 
body 1 

Suspension line spring constant, lb/ft 
Riser spring constant, lb/ft 
Length 

Length from confluence point to parachute CP 
Length of riser 

Length from SRB attach point to SRB CM 
Length of SRB 

ength from SRB CM to SRB CP, positive towards engine 
Length from confluence point to parachute CM 
Length of suspension lines 
Mass of body i 

Included mass of the parachute 
Apparent mass tensor of parachute 
Canopy mass 
Suspension line mass 

Moments about axis i of body 1 i = X Y, Z 
Moments about axis i of body 3 i = X, Y, Z 
Number of suspension lines, Normal force 
Angular Velocity about X-axis, body i 
Angular Velocity along Y-axis, body i 
Angular Velocity about Z-axis, body i 
Skirt Radius 

o 

Dynamic pressure at CP of body i, lb /ft 
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S oi 

2 

Nominal area, body i, ft 

u. 

X 

Linear velocity in X-direction, body i, ft/sec 

V. 

i 

Linear velocity ir Y-direction, body i, ft/sec 

w. 

1 

Linear velocity in Z-direction, body i, ft/sec 

X ., Y.. Z. 

Right-handed orthogonal axes of body fix d reference 
frame i 

X Ei’ Y Ei- Z Ei 

Earth fixed coordinates for body i, ft 

a. 

l 

Angle of attack, body i 


-1 V - 

Tan" , body i 

\lr 9., <t> . 

1 l i 

Euler angles, body i 

u 

Angular velocity vector 

C ) 

Dot notation for time derivative, 

C 

Damping coefficient, lb sec /ft 

P 

3 

Air density, slug/ft 


SUBSCRIPTS 

1 Parachute 

2 Riser 

3 Payload (SRB) 

o Nominal conditions 
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RECOVERY SYSTEM ANALYSIS 


THE APPROACH TO THE PROBLEM 

The technical approach is structured to assess the descent dynamcs and 
stability characteristics of a general parachute-riser payload combination. 
The advantage to a gene 1 case study is the ability to study a wider range 
of possible configurations with a minimum number of simplifying assump- 
tions. Three primary tasks describe the approach taken in the analysis of 
the problem. 

• The parachute riser-payload configuration was arranged 
and then said to be nominal according to specifications pro- 
vided by the contracting agency and particular requirements 
of the descent conditions. A mathematical model incorpor- 
ating elastic risers and suspension lines, three bodies each 
with six degrees of freedom, and a non steady air mass was 
developed. A complete softwa* c package was written to per- 
form the nonlinear simulation. 

• Using the nonlinear software package, simulations of the 
nonlinear dynamics of the parachute-riser-payload were 
made for a variety of initial conditions both with and without 
the influence of the nonsteady air mass and the elasticity of 
the suspension lines and riser. Particular attention was paid 
to equilibrium trajectories and to tl e occurrence of limit 
cycle responses. 

• Using numerical techniques, linearization of the state equa- 
tions of motion was accomplished. The stability of the 
system to disturbances was then assessed using the Root 
Locus technique. Using the same linearization technique, 
stability analysis as a function of certain parameters can 
be assessed. 

While the state of the art of parachute recovery of la: a e paylcads extrap- 
olates to a successful recovery of the space shuttle solid rocket booster 
(SRB), the magnitude of the SRB recovery problem is at least three times 
the size of any previous successful recovery. The large suspended load 
(approximately 150,000 lb), the size of the parachutes (3-130 ft Conical 
Ribbon) and the overall length of the system (about 400 ft) demand highly 
sophistocated math modeling and simulation if accurate stability conclusions 
are to be rightfully drawn. The ti_ hnical objective of his study then is to, 
as accurately as possible, analyze the descent dynamics, predict stability 
characteristics, and reduce the cost of the recovery by providing a better 
starting point for full scale testing and evaluation. 
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GENERAL RECOVERY SEQUENCE 


A schematic representation of the space shuttle SRB recovery is shown as 
Figure 1. The recovery process begins with the disengagement of the space 
shuttle main body and the SRB by explosive charges. The SRB then continues 
on a ballistic trajectory modified by its own aerodynamics through the 
apogee of nearly 200, 000 ft, descending to approximately 20, 000 ft, at which 
point the drogue parachute (48 ft Conical Ribbon) is deployed. Stabilization 
through the next 6000 ft of the descent provides sufficient conditions for the 
deployment in reefed stages of a three-parachute cluster. The cluster of 
130-ft conical ribbon parachutes is fully deployed and fully inflated at an 
altitude of approximately 6000 ft. A steady descent concludes with water 
impact. 

The analysis of the descent dynamics is made during the final 6000 ft, during 
which the motion of the system is effected by a potentially non-steady air 
mass perturbed by gusts. The analysis begins at full inflation of the cluster 
and ends at water impact. 

The recovery system components, the drogue parachute, the main para- 
chutes, and the SRB were chosen to meet the requirements established by 
the contracting agency. The drogue was chosen as a 48-ft, 20-cleg conical 
ribbon parachute (Ref. 1). A cluster of 3-130 ft, 20-deg conical ribbon 
parachutes provides the required 80 fps descent rate during the final 4000 ft 
(Ref. 2). The dimensions and mass of the SRB have continually changed 
during this study. The dimensions and mass used, however, are repre- 
sentative and provide an adequate model of the final configuration. The SRB / 
Drogue combination is shown in Figure 2, and the SRB/Main is shown in 
Figure 3. 


SIMPLIFYING ASSUMPTIONS 


Several simplifying assumptions are employed which reduce the computa- 
tional magnitude without compromising the general nature of the problem. 
Others are made to improve the math models to the extent that the state of 
the art allows. 

• The Parachute is assumed to be axisymmetric and to have 
a fixed-shape canopy with elastic suspension lines. 

• The riser connecting the parachute and payload is elastic 
and transmits only axial forces to the attach points on the 
SRB and parachute axes of symmetry. 

• The SRB is a rigid, axisymmetric body. 


6 




7 


Figure 1. Nominal SRB Recovery Sequence 


Drogue 


D q = 48 ft 

D x = 0.77D 0 = 37 ft 

1 $ = 2.0 D 0 = 96 ft 

L x = 103 ft 

L CM = 102 ft 

L CP = 0 163 0 o 

L $ - 96 ft 

SRB 

0 3 = 11.8 ft 

L 3 = 81 ft 

L 3T = 157 ft 

L a = 0.0 ft 


System 


t 

"2 

L 


Dq = 48 ft 
231 ft 



Figure 2. SRB/ Drogue Baseline Configuration 



Main 

Dq - 130 ft 

Lj = 292 ft 

L cp = 0.163 0 0 

Lcm = 290 h 

D x = 0 . 75 D 0 = 94 ft 

SRB 

D 3 = 11.8 ft 

L 3T = 145 ft 

L 3 = 75 ft 

l 4 = Oft 

System 

L 2 = 67 ft 



Figure 3. SRB/Main Baseline Configuration 
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• The aerodynamic centers of pressure are constrained to 
remain on the axes of symmetry of the SRB and the para- 
chute but do not necessarily coincide with the centers of 
mass of those bodies. 

• The energy modification of the air mass caused by the 
movement of the parachute through it is represented by 
tensors of apparent mass and apparent moments of inertia 
and not considered for the SRB motion. 

• The separation distance between the SRB and the main para 
chutes is large enough to neglect forebody wake effects. 

• The non steady air mass is represented by a wind velocity 
field and a gust velocity field perturbation. 

• A flat earth is used for trajectory calculations. 


SYSTEM MODELING 


The mathematical modeling of the primary subsystems, the parachute, the 
riser, and the solid rocket booster is described in this section as used in 
the development of an analytical nonlinear simulation programming system. 
Modeling of the elastic elements and the non steady air mass is also 
described. 

The equations of motion of the three body system are written relative to a 
liat earth. The forces and moments on the parachute and SRB result from 
aerodynamics and gravity. The application of the aerodynamics into the 
equations of motion is discussed. 

Finally, in this section the techniques used to linearize the nonlinear motion 
and to perform a stability analysis are outlined. 


DEVELOPMENT OF A NONLINEAR DYNAMICAL MODEL OF 
THE PARACHUTE/ RISER /PAYLOAD SYSTEM 


The parachute /riser /payload system is modeled as a three-body, six- 
degree-of-freedom-each problem. Since the parachute and SRB are con- 
nected by the riser, the constrained system finally reduces to a 15-degree- 
of-freedom problem. 
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The differential equations of motion. — The equations of motion are 
developed in general terms first with no elasticity and a steady air mass. 
The effects of the inclusion of the elastic suspension lines is then included. 
The basic math model is adapted from Reference 3. 

Reference frames: The reference frames and their initial orientation 
are shown in Figure 4. 

Four right handed orthogonal reference frames are needed to specify the 
motions of the parachute (System 1), the riser (System 2), and the payload 
(System 3). 

Earth fixed frame: Origin O e is fixed on an assumed flat earth 
directly below the initial position of the SRB Center of Mass. Z g is direct- 
downward, Xj, is horizontal on the flat earth aligned in the vertical earth 
plane containing the initial SRB Center of Mass velocity vector, and Y £ is 
cross range to the right. 

Body-fixed, moving frames 1, 2, and 3: The origins of the parachute 
and payload (SRB), body-fixed reference frames are at the respective 
centers of mass, Oj and Og. axes are aligned with the axes of symmetry 

with Zj directed toward the parachute confluence point, Zg directed from 

the parachute confluence point to the SRB attach point, and Zg directed 

toward the engine end of the SRB. axes are aligned initially parallel to 

the vertical earth plane containing the payload center of mass initial 
velocity vector. 

Euler angles: The Euler angles 0., 0., ^ describe the orientation of 

the body-fixed reference frames with respect to the earth fixed inertial 

frame. The ordered rotations are A . about Z. followed by 0. about Y. and 

i l J l l 

then about as illustrated in Figure 5. 

The direction cosine matrix [B-*] transforms a vector in earth fixed refer- 
ence frame to the jth body fixed reference frame in the following manner: 

V. = [B^] V E (1A) 

Conversely, by premultiplying by [B**] - * 

V E = [B j]' 1 V. (IB) 
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Figure 4. Reference Frame Definition and Orientation for 
a 3 -Body Parachute Riser Payload System 


12 




The direction cosine matrix [B**] is such that its transpose is the same as 
its inverse; i. e. , 

[B j ] T = [B ^]" 1 


Hence, Equation (IB) can also be written 


= [B j ] T V. 


(1C) 


In terms of the Euler angles and the sequence f , 0, 0 , [B-*] is as follows: 


[B j ] = 


cos cos iIn 


cos 0 . sin it. 
J y 3 


sin 0. cos \lt. sin 0 . sin 0 . sin 0 . sin i|t. 


3 


'J 


3 


3 


3 


3 


- cos 0 . sin i|r. 
3 3 


+ COS 0. COS lit. 

3 y 3 


cos 0 . sin 0 . cos i|t. cos 0 . sin 0 . sin i|t 


3 


3 


3 


3 


3 


-sin 0 . 

3 

sin 0 . cos 0 . 
3 3 


COS 0. cos 0. 
3 3 


+sin 0 . sin i|t^ 


-sir cos ijt. 

J J 


Its elements are written B^ where i is the row number and k is the 
column number. 

The Euler angle rates are given by 

i|r. = (Q. sin 0 . + R. cos 0 .) sec 0 


3 


3 


3 


J 


J 


J 


0. = Q. cos 0. - R. sin 0 . 
3 3 3 3 3 


(3) 


0 . = P. + (Q. sin 0 . + R. cos 0 .) tan 0. 
3 3 3 3 3 3 3 


The indices j = 1 , 2 , 3 correspond to the parachute, riser, and payload, 
respectively. 

The dynamics of motion: force and moment equations. -- 

The parachute: The equations of motion for the parachute are divided 
into force and moment equations about the center of mass. 
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The force equations are written 


F 1 + mjCB 1 ] g + [B 1 ] [B 2 ] T F, 





is the apparent mass tensor resulting from the air mass accelerations 
produced by the parachute motion. 
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Equation (4), when written in matrix form, becomes 




Equation (4), in scalar form, is 

F 1X + m l B 13 g + B S1 F 2 ' (m ! + m lAX ) <B 1 + W 1 ^1 ■ V 1 R l* 

F 1Y + m l B 23 g + B S2 F 2 = (m l + m lAY> < F 1 + U 1 R 1 * W 1 P l> < 4C > 


F 1Z + “1 B 33* + B S3 F 2 * <“l + “1AZ^ < W 1 + V 1 P 1 ' U 1 «1> 



1 2 T 

are the elements of the third column of the matrix [B ] [B ] 
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The aerodynamic forces are given by 


F 1X ■ 

C N1 * q l S ol )cos ^1 

F 1Y = 

C N1 * q l S ol )sin ^1 

F 1Z = 

- C T1 (q l ^l) 

-1 V 1 

where = tan (*jj )• 


The moment equations about the parachute body axes 
Mass may be written 

*i 

1 -♦ -4 

= hj + (^ x hj 


(5A) 


where = 


"1 

Ql 

R, 


the total angular velocity vector of body 1 and hj is the 


angular momentum vector of body 1 which can be written 



! XX1 0 0 


_P 1 ' 


P 1 r xxi 

hj = [I] u> 1 = 

0 I YY1 0 


Qi 

= 

Q 1 ! YY1 


0 0 *ZZ1 

— 




_ R 1 ! ZZ1_ 


The apparent moments of inertia resulting from the air mass accelerations 
generated by the parachute rotational motions may be written assuming 
principal axes 


I XXA1 ° 


= 


o 

o 


i 


YYAl 
0 I 


ZZA1 
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A combined moment of inertia matrix may be calculated, using the parallel 
axis theorem, and is written 


[1*1 


! XX1 

0 


0 

I* 

l Y* 1 


0 


0 


0 

0 

I* 

L ZZ1 


Hence, the moment equation may be written 


Mix 


• * 
*XX1 

0 

0 


"V 

m iy 

= 

0 

* 

T YY1 

0 



M 1Z 


0 

0 

* 

*ZZ1 


R 1 

_ 




- 


- 


0 




Q 1 


0 


-P 


1 



X XX1 
. 0 

0 


0 0 


P 1 

r YYl 0 



0 ! ZZ1 


R 1 


In scalar form. Equation (5B) becomes 


* 

M = I 

m lX l XXl 

P 1 + (I ZZ1 


* 

M = I 

iVA l Y l YYl 

Q 1 + (I XX1 

- ♦ 

* , ★ 

M 1Z = r ZZl 

R 1 + (I YY1 


- I YY i ) Q i R i 

- *ZZ1 > R 1 P 1 


- I 


XXI 


) p i Qj 


(5B) 


(5C) 


The moments acting about the CM location due to the external forces are 
in vector notation: 


M. ■ F 1 xL+[[B 1 ][B 2 ] T F,1xL r 
11 1 


(6A) 
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where 



In scalar form. Equation (6/ becomes: 


M 1X = F 1Y (L 1' L C M ) " F 2 B S2 L C m 



(6B) 


(6C) 
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There are no external forces acting off the axis of symmetry, hence 


The moment equations can be written using moment coefficients for the 
contribution to the total external moment due to aerodynamic forces as in 
Equation (74). 


The Payload (SRB) — The equations of motion for the payload are written 
along the same lines as those for the parachute, with the exception that the 
apparent mass and moment of inertia effects are not included. 


The force equations are written 

F 3 + m 3 [B 3 ] g - [B 3 ] [B 2 ] T F 2 = m 3 (C 3 +u 3 x Cg) 

where 


s = 

U 3 

V 3 

and 

u 3 = 

P 3 

q 3 


W 3 



R 3 


Equation (7), in scalar form, becomes 


F 3X + m 3 B 13 g ' B <m F *> = + Qq “ V, R,) 


S4 2 


3 w 3 "3 ^3 3 3 


F 3y + m 3 B 23 g - B qR F 9 = m^ (V a + U, R q - W, P ? ) 


'S5 2 


3 3 3 3 3 3 


F 3Z + m 3 B 33 g " B S6 F 2 ~ m 3 * W 3 f V 3 P 3 “ U 3 ^3* 


S4 

where ( B ct - ) are the elements of the third column of the matrix 

oO [ O n rp 

operation [B ] [B ] 


B 


S6J 


(7A) 


(7B) 
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The aerodynamic forces are given by 


P 3X = C N3 q 3 S o3 COS ^3 
F 3Y = C N3 q 3 S o3 sin ^3 
F 3Z = " C T3 q 3 S o3 

where 

-1 V 3 

^3 = tan < uf> 


The moment equations for the SRB are written 


M 

M 

M 


3X 

3Y 

3Z 


r XX3 P 3 + (I ZZ3 ' I YY3 ) Q 3 R 3 
r YY3 ^3 + (I XX3 " I ZZ3 ) R 3 P 3 
r ZZ3 R 3 + (I YY3 “ hcXS 0 P 3 Q 3 


(8A) 


The moments acting about the CM location due to the external forces are 

M 3X = " F 3Y L 4 ' F 2 B S5 L 3 
M 3Y = F 3X L 4 + F 2 B S4 L 3 
M 3Z = 0 

where L3 is the length from the SRB's attach point to its center of mass 
and L4 is the length f r om the center of mass to the center of pressure of 
the SRB* L4 is positive when measured from the center of mass in the direc- 
tion of +Zg. 

The moment equations can be written using moment coefficients for the 
contribution to the total external moment due to aerodynamic forces as in 
Equation (75). 


The Kinematics of Motion: The Riser Constraint 


The riser, assumed for the time being to be of fixed length, provides a 
convenient method of interconnecting the equations of motion of the parachute 
and the payload. Consider the linear velocities at each end of the riser. 
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At the confluence point of the parachute suspension lines: 



At the attach point on the payload: 



Subtracting Equation (9) from Equation (10), the linear velocities in the 
riser coordinate systems are eliminated: 



Differentiation of Equation (11) yields equations for Qg, P 2 and W^: 
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•dW 



' U l‘ 


f Q i L c 1 

1 C M 


1 T 

) - [B 1 ] < 


{t i 


^1 L C M ^ 
M 


V 1 

+ 

- p i B c 

1 M 



V 1 


- P 1 L C_. 

JVa j 


W 1 

1 

1 

0 




w 


0 

k 

L X J 

L 

- 

* 


k 1 

u - 


b 


( 12 ) 


The third scalar equation of Equation (12) gives an expression fo 1 W- as 
follows : 

^1 = TT* ”[®13 Q 2 L 2 + B 23 P 2 L 2 “ B 13 ^2 L 2 + B 23 P 2 T '2 
B 33 L 

+ B 13 (U 3 “ Q 3 L 3^ + B 23 (V 3 + P 3 L 3> + B 33 W 3 


+ B 13 ({j 3 - S L 3> + B 23 ( ^ 3 + P 3 L 3> + B 33 ™3 


(13) 


' B l 3 (U 1 + «l L C m > - B 23 (V 1 ' P 1 L C m > - B 33 W 1 
B 13 (U 1 + ^1 L C m ) " B 23 ( ^1 ' P 1 


Expressions for b|„. B 23' an< * B 33 ^ rom Equations (2), (3), and their deri- 
vatives are as follows: 


b] 


13 


" Q j B 33 + R j B 23 


Hi = P j ®33 ' R j B) 13 


(14) 


Hi • - p j Hi + B is 

Substitution from Equation (14) into Equation (13) yields expressions for Wj 
free of derivatives of B^ . . 

ik 

• * 

Similarly, we can obtain equations for Q 3 and P2 from the first and second 
scalar equations of Equation (12), 
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The riser force F 2 can be obtained from the third equation of Equation 
(4C) in terms of Wj as given by Equation (13). 

F 2 = B" f,m l +m l AZ ) [ V V 1 P r U l Q l 1 - F 12 - m l B 33 S 1 (15) 

The riser force, of course, is directed along the Z axis of the riser refer- 
ence frame. 


System State Differential Equations for the Non-Elastic, 

Steady Airmass Case 

Equations (3) to (8), (12), and (13) can be written in the following form: 


U. 


1 


m j + rrij 


f F lX +m l B 13 * + F 2 B S1^ ' W 1 Q 1 + V 1 R 1 < 16) 


AX 


V, 


m j + m j 


{F 1Y f mi B^ 3 g + F 2 B S2 1 - UjRj + WjPj (17) 


AY 


W. 


B 


I" ' L 2 f Q 2«2 B 33 - R 2 B 23> + P 2 ( P 2 B 33 ' R 2 B 13> 
33 


-Bj 3 Q 2 + B 23 P 2 )] + (R,B*, - Q,B^) <U, - Q, L,) 


3 23 

,3 


c 3 33 
>3 


c 3 3 


+B 13 * B 3 " < ^3 L 3 ) + ^ P 3 B 33 “ R 3 B 13 ) ' V 3 + P 3 L 3* 

+B 23 <^3 + P 3 L 3 ) + «3 B 13 ' P 3 B 23 ) < W 3> + B 33 < W 3> 

-( R l b' 23 - Qj B l 3 ) (Uj + Q l L c ^) - Bj 3 (&! + Qj L^) (18) 


■< P 1 B ls - R 1 B } 3 > (V 1 - P 1 L C m > - B 23 <*1 - P i 


-(Qi b } 3 - Pl b 1 23 ) (w x n 
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*3 ■ ^f F 3X + m 3 B ? 3 e- F 2 B S4l- W 3 Q 3 + V 3 R 3 
V 3 * nTT ' F 3Y + m 3 B 23 * ' F 2 B S 5 1 ' U 3 R S + W 3 p 3 


w, 


Q, 


R. 


Q, 


R. 


HTJ t F 3Z + m 3 B 33 « * F 2 Bse 1 ‘ V 3 P 3 + U 3 «3 


fxx7 f M lX " (T ZZ1 " r YYl ) Q 1 R ll 

I* Y1 t M lY " (I XX1 “ ! ZZ1 ) R 1 

I* Z1 f M lZ " (I YY1 " *XX1 ) P 1 ^ 

^XX3 ^ Msx " <lzz3 ~ I YY3 ) Q 3 R 3^ 


[M„v “ 


r YY3 3Y XX3 ZZ3' *3 3 


I ZZ3 f M 3Z " (I YY3 ‘ I XX3 ) P 3 Q 3^ 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 

(25) 

(26) 
(27) 


0. 

J 


(Q. sin <£. + R. cos <J>.) sec 0., i = 1, 2, 3 
J J J J J J 


Q. cos <t>. - R. sin <fl., j = 1, 2, 3 
J J J 3 


(28)-(30) 

(31)-(33) 


<$ . = P. + (Q. sin <A. + R. cos <j.) tan 0., j = 1, 2, 3 

J J J J J J 3 


(34)-(36) 


B 


1 («2 <Q 2 B n - R 2 B 22> + P 2< P 2 B 32 ‘ R 2 B ? 2 > ‘ «2 ®12 


22 


+ -q C(R 3 B 22 - % B 32> (U 3 ' «3 L 3> + B 12 <U 3 * «3 L 3> 
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+ (P 3 B 32 ' R 3 B 12> (V 3 + P 3 L 3> + B 22 <V 3 + P 3 L 3> 


+ <Q 3 B 12 - P 3 B 22> W 3 + B 32 <Wj> 


( 37 ) 


- < R 1 B 22 - «1 B 32> (U 1 + «1 L C m » - B l 2 <4 + «1 L C m » 

- <P 1 B 32 ‘ R 1 B l 2 » (V 1 ' P 1 L C > - B 22 1 - P 1 L C> 

M M 

- <«1 B 12 ' P 1 B 22> W 1 - B 32 

^2 * + 3“ f®2 ( ®2 B 31 " R 2 B 21* + P 2 <P 2 ®31 ' R 2 B ll’ + P 2 B 21 
B 11 

+ rr 2 t <R 3 B 21 - % B 31> <U 3 - «3 L 3> + B U <4 ‘ «3 L 3> 

+ < P 3 B 31 - R 3 B n> (V 3 + P 3 L 3» + B 21 ,V 3 + P 3 L 3> 


+ ,Q 3 B ?1- P 3 B 21 ,W 3 + B 31 ( ' i, 3 ) 


(38) 


- < R l ®21 - Q ’ B Jl> (U 1 + Q 1 L C „> - B ll <«1 + «1 L C, » 


M 


'M 


- (p i b 5i - R 1 B n> (v i - p i L C M > - B 2i <4 - p i L C M > 
B n ' p j B 2i* w i ' B 3i 


R 2 * 0 


(39) 


In the set of Equations (16) through (39), there are no terms involving the 
riser linear velocity components Ug, Vg, Wg or their derivatives. 

The velocity of the SRB Center of Mass relative to the earth can be deter- 
mined from Equations (19), (20), (21), the direction cosine matrix [B^] 
given by Equation (2), and Equation (1C). 
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By integrating Equations (19), (20), and (21) and applying Equation (1C), 
the linear velocity components of the SRB Center of Mass in the Earth 
fixed reference frame will be: 


E3 = U 3 B U + V 3 B 21 + W 3 B 31 

(40) 

E3 = U 3 B 12 + V 3 B 22 + W 3 B 32 

(41) 

E3 = U 3 B ?3 + V 3 B 23 + W 3 B 33 

(42) 


Elastic Models 

The entire parachute is made of elastic material and is subject to deforma- 
tion under load. The riser too elongates when loaded. To account for any 
additional dynamics caused by the continuous dynamic flexing of the para- 
chute-riser system, two elastic models are employed. 

The first is a canopy shape model which depends on the suspended load, 
the canopy pressure distribution, the inflation condition, and the construc- 
tion of the canopy (Ref. 4). Its use is independent of the simulation pro- 
gram but its result is an input to the simulation program. 

The second is a damped spring mass model of the suspension lines and 
riser. The application of this model is dynamic in the simulation program. 

Program CANO -- An elastic canopy shape analysis is done by Program 
CANO (Ref. 4). 

For an assumed pressure distribution and an initial gore geometry specifi- 
cation (which assumes the canopy to be made up of discrete horizontal and 
radial elements) and a specified suspension line length and riser load the 
program solves for the equilibrium shape and loads of the discrete mem- 
bers. 

The method assumes an elastic deformable frame (the canopy) under a 
specific load (the pressure distribution) to determine the loaded (equili- 
brium) shape. The pressure distribution is nondimens ionalized by the 
length along the canopy surface. The load elongation curves are set for 
types of materials and are generalized as percentage of breaking strength 
and unit strain. Thus, only the type of materials and geometry of the gore 
need be specified. 

For specific loading conditions such as reefed, fully open, overinflation 
lines, etc., the program iterates across the canopy surface, adjusting the 


27 



breaking strength of each member to equal the calculated load. This in 
turn adjusts the weight and then the equilibrium shape of the canopy. For 
the new breaking strengths of the members, a new equilibrium shape for 
the canopy and new loads for each element are calculated. The loads are 
compared to the most currently defined breaking strength, and when ail 
elements have breaking strengths within a range of zero to five percent 
more than the calculated load, the parachute is said to be optimized. 

Using the assumption that for a particular material type (e. g. , web, tape, 
or cord) the weight of a material is proportional to its breaking strength, 
the optimized weight of ihe radial and horizontal members and the suspen- 
sion lines are calculated, and thus the weight of the total canopy is deter- 
mined. From an input table of available materials characterized by break- 
ing strength and type, which implies a load-strain characteristic and a 
parametric weight, materials are chosen that are the lightest available 
which meet the strength requirements for the calculated loads in individual 
elements. The "buildable" parachute weights are calculated and compared 
as non-optimum factors to the optimized parachute weights. 

The program CANO can be applied to consecutive steps in the process of 
deployment. It can be used to calculate the optimum parachute to meet up 
to 21 loading conditions which are combinations of partial inflation, reefed 
skirt, overinflation lines, and fully open. Thus, an accurate estimate of 
canopy weight can be made for a particular set of loading and inflation 
conditions. 

The canopy profile generated by CANO for a fully inflated 130 ft conical 
ribbon parachute with a 200, 000 lb suspended weight is shown in Figure 6. 

Elastic Suspension Lines and Risers -- The suspension lines and risers 
generally used in parachute construction are quite elastic. The additional 
dynamics introduced by their elastic characteristics are to be included in 
the general equations of motion describing the parachute /riser /payload 
descent. 

The geometry of the parachute and riser is shown in Figure 7. The elastic 
elements are the suspension lines (length L s ) and the riser (length L 2 ). 
Elongation of the suspension lines results in a change in che suspension 
line angle and hence the suspension line moments of inertia. There is also 
a change in the location of the center of mass of the parachute and a result- 
ing change in the total moments of inertia of the parachute. 

Two key assumptions are made: 

• The canopy is fixed in shape and thus the skirt radius R 

is constant. ° 

• The angle between the parachute axis of symmetry and the riser 
is always small. The riser force then is transmitted to the 
confluence point along the parachute axis of symmetry and thus 
the parachute remains axially symmetric. That is, the con- 
fluence point remains on the axis of symmetry and the suspen- 
sion line cone remains right circular with variations in height 
only. 
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The elastic elements are modeled as damped linearly elastic springs. The 
damping coefficient is taken as a representative value for dacron material 
(Ref. 5). Thus, 


c = 0.05 lb sec /ft (43) 

The spring constants are determined as functions of the unstretched length, 
the elongation at break, the suspended load, and a safety factor of 3. Thus, 
for the riser. 



3M 3 g 
1. 2 L t 


lb/ ft 


(44) 


where Mgg is the suspended weight 

is the unstretched riser length 
o 

1. 2 represents 20% elongation at break 
and 3 is a safety factor 

For the suspension lines the suspended load is the load carried by each 
line so that 


K 


LS 


3M 3 g 

2L So^ 


lb /ft 


(45) 


where Lg 0 is the unstretched suspension line length and N is the number of 
suspension lines. 

We can now model the dynamic length of the elastic elements for the sus- 
pension lines 


= u 


So 


L <V ttn 'l ( S2) 

N K l cos y ’ y xan 


for the riser 


J R = L Ro + ' 


(F 9 - C l 2 > 

Kr, 


(46) 


(47) 


♦ » 

In Equation (46) Lqm is used in place of Ls. Differentiation of Equation 
(53) (yet to come) with respect to time validates this substitution provided 
that the suspension line angle is small and the included mass much greater 
than the canopy or suspension line mass. 
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The rate of change of lengths of the center of mass location and the riser 
are calculated using the central difference in average length divided by the 
change in time. Thus, 


dL 2 


L 2 - l 2 | 

* t- At -* t It- 

dt 

t-At 


2 At 
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4| t is the averaged riser length during the interval 
is the time derivative of at time t 
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(48) 


Parachute Center of Mass Location -- The canopy is modeled as a semi- 
oblate spheroid whose height is 32. 5% of the nominal diameter, D 0 and 
whose radius is 36% of the nominal diameter. 


The canopy volume is 


V = | tt (0.325 D ) (0.36 D ) 2 . 
C 3 O O 

The included air mass is given by 

m I V c» 


(49) 


(50) 


where p is the air density. 
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The center of pressure and the center of mass of the canopy are both 
assumed to be located at the centroid of the canopy volume. 

The suspension lines are modeled as thin rods having uniform mass dis- 
tribution as shown in Figure 8. 



Figure 8. Suspe . i '-e Geometry 
The suspension line angle is given by 


y 


tan 



The center of mass of all the suspension lines is 


AJLCM _±S 


cos y 


( 51 ) 


( 52 ) 


The cai ter of mass location for the entire canopy then is given by the fol- 
lowing relation 


L 



ALCM x m 


+ L, x m ■) L, x m. 

1 c 1 I 

m, + m + m r 
1 c l 


( 53 ) 
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Parachute moments of inertia-. Canopy moments of inertia written about 
the point "o" in Figure 6 in the plane of the skirt are 


>ZZ = T M c ( 0 - 36 D o )2 


1 = I = 4 M [(0. 325 D ) 2 + (0. 36 D ) 2 ] 

AA YY 3 C O O 


(54) 


Suspension line moments of inertia about the suspension line jone center of 
mass location are 


l ZZ 


M ^ L S 

12 


*xx 1 


. 2 
sin y 




(55) 


YY 


12 


2 

cos“y 


The apparent moments of inertia of the canopy written about the total para- 
chute center of mass location are according to Reference 6: 

I„ 7 = 0.063 p(R) 5 

zz Aj ° 


‘xx ' >YY Ai ■ 0.042 p(R o ) 5 + M 1 a ( Li -L Cm ) 2 


(56) 


The total Parachute moments of inertia about the total parachute center of 
mass location then are given by 
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Additions to the nonlinear differential equations of motion, -- Tin; in- 
clusion of elasticity adds several terms to the differential equations of 
motion. 


When writing the constraint equations which allow coupling of the motions 
of the parachute and SRB, the velocity of the confluence point relative to 
the center of mass location is amended to read 
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(58) 


The velocity of the end of the 


riser at the payload attach point is rewritten 
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~«2 L 2 
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+ 
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4 


(59) 


Thus, with the addition of elasticity, the final constraint equation, corre- 
sponding to Equation (11) is: 
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(60) 
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Differentiating the constraint equation above results in 
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From Equation (61) we can obtain expressions for Wj , Pg, Qg* 

Equations (16) to (39) describe the nonelastic differential equations of 
motion. Equations for Wj, P 2 , and Q 2 C(18), (37), and (38)] are re- 
written here to incorporate the changes due to the inclusion of elasticity. 
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DEVELOPMENT OF AN AERODYNAMIC FORCE AND 
MOMENT SYSTEM IN AN UNSTEADY AIRMASS 


The application of aerodynamic forces and moments in the simulation pro- 
gram is described in this section along with the models describing the non- 
steady air mass and their effect on the aerodynamic forces and moments. 

Least squares polynomial curve fits to the aerodynamic coefficients data as 
given in References 9 and 10 provide a convenient method of representation 
of the normal, tangential, and moment coefficients as functions of the angle 
of attack for the parachute and SRB. The parachute force and moment sys- 
tem is shown in Figure 9 followed by the normal and tangent force and 
moment coefficients curve fits in Figures 10, 11, and 12, respectively. 

The SRB force and moment system is illustrated in Figure 13 and the SRB 
normal and tangent force and moment coefficients curve fits in Figures 14, 
15, and 16, respectively. 

Wind and gust models. -- To determine water entry characteristics of 
the SRB, the effects of winds and gusts near the surface of the earth on the 
attitude of the descending SRB must be accounted for. 

Wind and gust models to provide inputs to the recovery simulation as re- 
quired by the contracting agency are described as adapted from Reference 7 
and Reference 8, respectively. 

Wind model: The recovery analy of the space shuttle Solid Rocket 
Booster <SRB) requires steady-state winds to be defined in the layer of air 
between sea level and 3281 ft (1 KM). The following is the recommended 
5% risk steady-state wind profile of Reference 7. 

V wind (h) ' V wind (h Q ) -str” 0 ih * 582 ft 

V wi„d (h) 1 V wind (h > 582 Sh * 3281 « 

o 

(65) 

V . , ,, , = 69 fps 

wind (h ) r 

o 

h = 3281 ft 
o 

P = 0.21 

The steady- state wind profile is shown in Figure 15. 


39 





l/> 


1 


o 


<3 


cm in 
r** o 

rH r* 

cm co 


O rH 

x x 



i 

UJ 

O 

co rj 
co <n 
o £ 

• CO 


o 


UJ 

CO 

CM co 2 

x x oc 


o 


41 


Figure 10. Polynomial Curve Fit to Normal Force Cc 
cients for a 20° Conical Ribbon Parachute 
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Figure 12. Polynomial Curve Fit to Moment Coefficients for a 70 
Conical Ribbon Parachute 






Figure 14. Polynomial Curve Fit to Normal Force Coefficients 
for the SRB 
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Figure 16. Polynomial Curve Fit to the Mom it Coefficients of the 



Gust model: Associated with the steady-state wind profile (air mass 
velocity field) is a discrete gust environment. The gust amplitude repre- 
sents a step velocity change in the air mass velocity field. The maximum 
gust amplitude envelope associated with the 5% risk steady- state wind pro- 
file as recommended by Reference 8 is as follows: 

V = 19. 7 fps 0 * h * 980 ft 

gust K 

v * = %Srr (h-980) + 19.7 980 S h s 3281 ft (66 (66) 

gust 2301 

V . = 29.7 fps h > 3281 ft 

f iSt 

The gust envelope is superimposed on the steady-state wind profile in Fig- 
ure 17. 

Relative velocity vector: The aerodynamic forces and moments are 
functions of the angle of attack, the altitude, the nominal area, a reference 
length for the moments, and the velocity vector of the center of mass with 
respect to the wind. 

The velocity field of the moving air mass can be written 

V = V . . + V . 
wg wind gust 


whe re 


V is the velocity field vector 
wg J 


V . , is the mean wind velocity field vector 

wind 17 


Vgust * s S us * vel <> cit y field vector 
and all are, in general, altitude dependent. 

The influence of the motion of the air on the body aerodynamics is accounted 
for by determining the velocity of the body with respect to the air to be used 
in developing the aerodynamic forces and moments. The relative motion of 
the center of pressure appears then as 

V a = C + wxL - Cb 1 3 V wg (67) 
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where 


C is the velocity of the center of mass with respect to the 
earth in body coordinate directions 

u is the angular velocity of the body 

I- is the vector from the body CM to the body CP 

Written in matrix form, for body i 



These are the velocity components used to determine the aerodynamic 
forces and moments. 

Elsewhere in the dynamical equations, inertial velocities are used. 
The angle of attack is given by 
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The side slip angle is defined as 


0. = Tan’ 


(70) 


X. 

1 


Aerodynamic Forces and Moments 

The aerodynamic forces acting on the parachute can be written in the body 
fixed axes directions. 


(71) 
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Similarly for the SRB the aerodynamic force in the body fixed directions are 

*3 


F o 

3 X 


F 3 

3 Y 
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C w Cos 0, 
N 3 


C N 3 Sin ^3 
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0 0 
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’^3 S o- 


(72) 


In general aerodynamic moments are written in terms of a moment coeffi- 
cient (Cm) and a reference length (MRP). The aerodynamic moments then 
about the X and Y body fixed axes whose origin is located at the MRP can be 
written. 


(73) 
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The moment reference point length (MRP) for parachutes is generally one 
nominal diameter ahead of the skirt plane. 

To write the aerodynamic moments about the body fixed axes system lo- 
cated at the body center of mass, new MRP lengths must be defined. 


The normal force is experimentally measured at the vent of the parachute. 
The height of the canopy plus the moment reference length is given by 

0. 325 D +D = 1. 325 D 
o o o 


The moment then is 


N(l. 325 D q ) , where N is the normal force 
The distance from the vent to the center of mass of the parachute is given 

by 


0. 325 D + ALCM - LCM 
o 

The functional form then of the aerodynamic moments acting on the para- 
chute written about the body fixed axes located at the parachute center of 
mass is 


M 1 


’ C Mj Si " *1 

M 

1 Y 


C Mj CoS ^1 



» - 


D (0. 325 D + ALCM - LCM) 
_ o o o 

* ' q l & o. 1.325 D o 
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D (0. 325 D + ALCM - LCM) 

q l S o. 1.3S5 D 

1 o 


The SRB aerodynamic moment coefficients are defined by 



N (L 3 + L 4 ) 


2 ’ V 53 ° 3 


where N is the normal force 
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The functional form of the aerodynamic moments acting on the SRB written 
about the body fixed axes located at the SRB center of mass due to aero- 
dynamic normal forces acting at the center of pressure is 


1 

s 

03 

* 1 


C M 3 Sin 3 3 


D 3 L 4 

< 4+l* 0 

M 3 


C M 3 Cos <*3 
— - 



L y J 



0 ,s 

“3% ( L 3+ L 4 ) j 


(75) 


LINEARIZATION OF THE EQUATIONS OF MOTION 

Application of the root locus stability analysis techniques to the solution of 
the SRB recovery problem requires a linearized system of equations of 
motion. One method of linearization is to choose a reference state, say 
vertical descent, and define small disturbances about this state. After 
linearizing the aerodynamic coefficients with respect to small changes in 
angle of attack and making appropriate substitutions, the linearized state 
is obtained by neglecting terms of order 2 and higher. This is a cumber- 
some task and the result is applicable only to the particular reference 
state originally chosen. 

A more general linearization method results from numerical techniques 
developed in Reference 11. 


Linearization Technique 


For a nonlinear system of equations implicit in time, the state can be re- 
presented as 

• » 

x = f (x, x) (76) 

where 


x = x (t) 
x = £(5(t)) 
t = time 


We want to linearize the vector nonlinear differential equations represented 
by Equation (76) at a particular point in time t Q . 


53 



The methodology is to calculate the nonlinear solution of x until t = t„ and 
then use the nonlinear solution at t 0 as the reference state about whicn the 
equations of motion are linearized. 

Let x be the nonlinear solution of Equation (76) at time t Q and 3E be the 

linearized solution at t : 

o 

x is known 

x is to be numerically derived 
more explicitly 


F^(x, x) 


g|- q (t, x) ..... f, (x. x) 

• • 

* • 

• * 

W* f n (X * m 55~ f n (X ' 50 

l n 


(77) 


Here 
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*1 


*1 

f 2 


x 2 

f 3 

X = 
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5E 

n 


n 


for a system of n equations. 

The matrix F x (x, x) represents the first partial derivatives of each state 
equation with respect to each state variable. The elements of F x (x, x) 
are determined by the central difference quotient 


a3C j 1 ^5 J7 8 ) 


where A*- is taken to be 1 percent of X.. 

J j 
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Further, let ? = *x - x be the disturbance vector about x. Differentiation 
yields 


- i, 


Rearranging terms 


? = x-x = f (x, x) - f (x, x) 


X = X - 5 


(79) 


? = f(x + e, x + ?) - f(x,x). 

The mean value theorem of differential calculus allows 


? = f(x + $, x + ?) ~ F-(x,x) ? 


(80) 


where 




and AXj is the disturbance of the element, x^, of x. 


Equation (80) can be solved using the matrix of partial derivatives (77). 

The solution, call it y, is linear and the desired linearized state is found 

x = x + y. (81) 


Eigenvalues 

th 

Manipulation of the coefficients matrix of Equation (80) results in an n 
degree characteristic polynomial whose n roots are the eigenvalues. 

Actually the solution to Equation (80) is not found because only the eigen- 
values are required. The matrix of system (80) is transformed to Upper 
Hessenberg form. Using a Q-R procedure with double iterations and a 
convergence check, the eigenvalues to Equation (80) are approximated. 

The eigenvalues are of the form 

<7 ± jW 

where a is the real part 

u is the damped frequency 

j isV-T 
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STABILITY ANALYSIS TECHNIQUE 


The Root locus technique plots the eigenvalues on a complex plane. The 
relative stability and transient performance of the system are directly re- 
lated to the position of the eigenvalues. The root locus plot provides a tool 
for investigating the effect of parametric variations on system response 
and stability. The sensitivity to adjustments of a particular parameter can 
be examined and a systematic procedure can be followed to move the root 
locus to a desired position on the complex plane corresponding to required 
stability and response characteristics. 


ANALYSIS OF THE SOLID ROCKET 
BOOSTER RECOVERY SYSTEM 


To determine an entry envelope of orientations of the SRB as functions of 
initial conditions, elasticity dynamics, and nonsteady air mass conditions, 
a wide variety of simulations were made on the nominal descent configura- 
tions from an altitude of 6000 ft to water impact after approximately 74 
seconds. 


NOMINAL BASELINE CONFIGURATIONS 

The drogue and main parachutes in combination with the SRB were illustrated 
in Figures 2 and 3, respectively. Their specific dimensions are listed 
in Table 1. 


SINGLE PARACHUTE EQUIVALENCE TO THE CLUSTER 

The cluster of parachutes is modeled by a single parachute having the physi- 
cal dimensions of one of the parachutes in the cluster but the mass, inertia, 
and drag area characteristics of the entire cluster. 

In program CHUTER, described in Appendix A, all of the parachute- 
related input data are for a single element of the cluster. The number of 
chutes in the cluster is also a data input. The conversion to the equivalent 
parachute is handled within the program. 


NOMINAL SYSTEMS RESPONSE TO DISTURBANCES 

Two principal modes of disturbance or initial conditions were used in exam- 
ining the nominal systems response to initial conditions. For analytical 
purposes, the disturbances are induced in only one plane and thus the motions 
are in one plane only. A "pendulum" disturbance in which the parachute. 
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TABLE 1 - RECOVERY SYSTEM PARAMETERS 


Drogue/ Main/SRB 

SRB (Equivalent) 


Parachute 



o 

Q 

48 ft 

130 ft 

S 

°1 

1810 ft 2 

39900 ft 2 

L s 

96 ft 

275 ft 

1 L 1 

100 ft 

310 ft 

M c 

11 slugs 

69. 9 slugs 

m l 

9 slugs 

81.6 slugs 

Riser 



L 2 

48 ft 

67 ft 

SRB 



D °3 

11.8 ft 

11. 8 ft 

% 

110. 0 ft 2 

110. 0 ft 2 

T 

“3 

81 ft 

75 ft 

L 3„, 

T 

157 ft 

145 ft 

m 3 

5000 slugs 

4750. 0 slugs 

>xx 3 

8. 36 x 10 6 

7.36 x 10 6 

'yy 3 

8. 36 x 10 6 

7. 36 x 10® 

'zz 3 

1.96 x 10 6 

1.72 x 10 6 










riser, and payload remain generally aligned while being tipped to some 
initial angle results in smaller angular excursions of the SRB with less 
damping in the transient phase of the response. 

A "scissors" disturbance is one where the parachute and riser are marked- 
ly misaligned with the SRB. Response to this initial condition results in 
larger SRB angular excursions but with higher damping in the transient 
phase. 

Several sets of each type of initial conditions were imposed on the SRB / 
Main parachute combination. To see the added effects of elasticity and 
wind, each set was first run without the elastic or nonsteady air mass op- 
tions. The same cases were then run with the addition of elasticity only 
and rerun again with the nonsteady air mass option only. 

For reference, a case with no initial disturbance was run without elasticity 
or wind, with wind only, and with elasticity only. 

The cases specifically illustrated are listed in Table 2. 

TABLE 2 - ILLUSTRATED NONLINEAR 
SIMULATION CASES 


Configuration 

Initial 

Displacement 

Type 

01 

(fleg) 

03 

(deg) 

Elastic 

Winds 

and 

Gust 

Nonlinear 

Response 

Figures 

Root 

Locus 

Figures 

SRB/Main 

Pendulum 

+20 

+20 

No 

No 

18, 19 

41, 42 

SRB/Main 

Scissors 

-20 

+20 

No 

No 

21, 22 

43, 44 

SRB /Drogue 

See Fig. No. 24 

--- 

--- 

No 

No 

25, 26 

— 

SRB/Main 

Vertical 

0 

0 

No 

Yes 

20, 29 

45, 46 

SRB/Main 

Pendulum 

-20 

-20 

No 

Yes 

30, 31 

— 

SRB/Main 

Pendulum 

+20 

+20 

No 

Yes 

33, 34 

47, 48 

SRB/Main 

Scissors 

+20 

-20 

No 

Yes 

35, 36 

... 

SRB/Main 

Pendulum 

+20 

20 

Yes 

No 

37, 38 

49. 50 

SRB/Main 

Scissors 

-20 

+20 

Yes 

No 

39, 40 

--- 


SRB /Main Parachute Response to Pendulum -Type 
Initial Displacements 

Pendulum -type initial disturbances of up to 30 degrees were imposed on 
the SRB/Main parachute descent configuration. The responses were similar 
in nature so that only the angular response for a +20 deg pendulum-type 
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disturbance is shown as Figure 18. In all pendulum-type initial distur- 
bances with a steady air mass the parachute angular orientation over shoots 
by approximately 45% and the SRB angular orientation over shoots by approx- 
imately 55%. The response is typified by the shorter period oscillations of 
the SRB as it follows the orientation of the parachute. The relative motions 
of the parachute and SRB quickly become 180 deg out of phase, and the SRB 
motion induces perturbations on the long period parachute response. The 
angle of attack time history is depicted in Figure 19 and the trajectory is 
shown in Figure 20. 


SRB/Main Parachute Response to Scissors-Type 
Initial Conditions 


Scissors-type initial conditions of up to 60 deg misalignment were imposed 
on the SRB/Main parachute descent configurations. The responses for a 
scissors-type displacement with no wind or elasti ty were similar so that 
the angular response for only one parachute initial angular displacement 
of -20 deg and a SRB initial angular disturbance of +20 deg is shown in 
(Figure 21). Scissors-type initial conditions produced responses typified 
by 180 deg out of phase oscillations of the parachute and SRB, with the 
parachute motion, again long period, driving the general motion of the SRIi 
and the SRB inducing small perturbations on the otherwise smooth para- 
chute response. As in the cases with pendulum displacements the parachute 
over shoots to approximately 55%. The greatly increased moments on the 
SRB cause overshoots of approximately 170%. The long-term result of a 
scissors displacement is larger SRB angular excursions through the entire 
descent. The angle of attack time history is shown in Figure 22 and the 
trajectory is shown as Figure 23. 


SRB/Drogue Response to an Assumed Deployment 
Condition 


The SRB, after leaving the space shuttle, is assumed to move along a tra- 
jectory with a large angle of attack near 90 deg. Additionally, the SRB 
may be spinning about an axis approximately parallel to the trajectory. 

The object of the drogue parachute is to stabilize the SRB; that is, reduce 
its angle of attack to sufficient conditions required for deployment of the 
main parachutes. If the SRB is spinning, the drogue parachute will also 
reduce the total angular velocity of the SRB. 

The SRB/Drogue combination is simulated at an altitude of 20000 ft descend 
ing vertically at a rate of 580 fps. Its initial angle of attack is taken to be 
80 deg and the SRB is assumed to be rotating at 40 deg/sec about the earth 
fixed Z axis. The drogue parachute, assumed to be previously deployed, 
is initially positioned at a 10 deg yaw angle. The initial conditions are 
illustrated in Figure 24. 
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Figure 18. Main Parachute - SRB Response to a 20-deg 

1 endulum Disturbance (X-SRB, 2-sec Intervals) 




Figure 19. Main Parachute - SRB Angle of Attack Response 
to a 20- deg Pendulum Disturbance (X-SRB, 
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Figure 20. Recovery System Trajectory, Pendulum Initial 
Conditions 




Figure 21. Meta Parachute - SRB Response to a Scieeors 
Displacement ( 9 1 « *20 deg, 9 3 * +20 deg 
Initially) (X-SRB, 2-aec Intervals) 





ire 22. Main Parachute - SRB Angle of Attack Response 
to Scissors Displacement (0 1 « -20 deg, 0 3 * 
+20 deg Initially) (X-SRB, 2 sec Intervals) , 
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Figure 23. Recovery System Trajectory, Sciesors 
Initial Conditions 




The coning angle is a combination of the Euler angles 9 and 9 

cone angle = cos" 1 (cos 9 cos 9) 

It is the angle between the vertical descent line and the axis of symmetry 
of the body. The reduction of die SRB cone angle by the action of the drogue 
parachute is shown in Figure 25. The angle of attack time histories of the 
parachute and SRB are shown in Figure 26. 


Additional Effects Due to a Steady Wind and Gusts 

The application of an air mass velocity profile (mean wind plus gusts) as 
shown in Figure 27 to the descending SRB/Main parachute configuration, 
which is previously undisturbed, causes a rapid increase in the downrange 
velocity of the entire system. Figure 28 shows the Euler angle, theta, time 
history of die parachute, and SRB whose initial conditions were vertical 
descent. For the same case Figure 29 shows the angle of attack time his- 
tory. The initial large positive angle of attack produces large normal 
aerodynamic forces on the parachute. The parachute swings to a large 
negative orientation angle. The SRB, with a shorter period, being driven 
by the motion of the parachute, again follows. The parachute angle of 
attack quickly reduces to small angles while the SRB with far less aero- 
dynamic pitch damping requires more time to stabilize and damp its angle 
of attack. 

Pendulum Initial Conditions — Since the parachute is the driving force in 
the motion of the recovery system, its orientation initially with respect to 
a nonsteady air mass dictates the system response. Figure 30 depicts the 
Euler angle theta, time history for the SRB/Main parachute recovery sys- 
tem tipped down wind at -20 deg. The SRB and parachute orientation angles 
respond quickly to gusts at 15 sec and 45 sec. The overall response in the 
nonsteady air mass is stable. The angle of attack time history for the 
down wind pendulum case is shown as Figure 31. The gust can easily be 
seen as large sudden changes in the angle of attack. The SRB angle of 
attack decreases near the ground as the air mass velocity field slows down 
in the boundary layer effect. 

A trajectory typical of all cases run with nonsteady air mass is shown in 
Figure 32. 

If the parachute and SRB in a pendulum displacement mode are tipped into 
the wind, the response, although similar to the pendulum displacement 
downwind, is more dramatic. The increased angular excursions for a case 
tipped +20 deg is seen in Figure 33. Similarly, while the characteristic 
shape of the angle of attack time history for pendulum initial conditions is 
evident, the increased amplitudes for the system tipped into the wind initially 
are evident in Figure 34. 
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Figure 26. Drogue Parachute /SRB Angle of Attack Response 
to Initial Conditions Illustrated in Figure 24 
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Figure 29. Main Parachute - SRB Angle of Attack Response 

to a Non Steady Air Mass (X-SRB, 2-sec Intervals) 



MAIN PARACHUTE. PAYLOAD 
ORIENTATION T I HE HISTORY 



Figure 30. Main Parachute - SUB Response to a Downwinc 
Pendulum Initial Condition In a Non Steady Air 
Mass (X-SRB, 2-sec Intervals) 
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Figure 31. Main Parachute - SRB Angle of Attack Response 
to a Downwind Pendulum Initial Condition in a 
Non Steady Air Mass (X-SRB, 2 -sec Intervals) 




Figure 32. Recovery System Trajectory in a Non Steady 
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Figure 33. Main Parachute - SRB Response to an Upwind 
Pendulum Initial Condition in a Non Steady Air 
Maas (X-SRB, 2- sec Intervals) 
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Figure 34. " T ’ in Parachute - SRB Angle of Attack Response 
to n Upwind Pendulum Initial Condition in a 
Non Steady Air Mass (X-SRB, 2-sec Intervals) 


Scissors Initial Conditions — As seen in Figure 21 the scissor mode initial 
conditions result in larger amplitude SRB oscillation. The application of 
an altitude-air mass velocity profile as in Figure 27 to scissors mode 
initial conditions of the SRB /Main parachute primarily causes a down wind 
drift approximately equal to the wind speed. During the initial transient 
period when the recovery system is accelerating down wind, large angular 
excursions of both the parachute and SRB are seen (Figure 35). The SRB 
ang’e of attack becomes quite large as seen in Figure 36. The stability of 
tl.i system is evident at 15 and 45 sec as seen in the angulai response 
(Figure 35) to gust inputs. 


Additional Effects Due to Elasticity 

The inclusion of the elastic suspension line model in the nonlinear simula- 
tion allows the geometry of the system to be dynamically variable. The 
change in suspension line lengths in particular changes the mass distribution 
of the parachute slightly; thus, through the change in moments of inertia a 
slight decrease in the period of th*~ parachute is seen. 

In the differentiated constraint [Equation (61)] which includes the elastic 
suspension system, the velocities and accelerations betweer the u.id points 
of the riser and the confluence point and center of mass location are re- 
quired. The elastic elements flex at several frequencies depending on the 
frequencies of the parachute, riser, and SRB oscillations. To calculate 
the velocities and accelerations required, a nume rical method was used to 
average the lengths over the high-frequency oscillations and then calculate 
the rates based on a frequency approximately one-half of the SRB natural 
frequency. This frequency was chosen since the riser force peaks at each 
local maximum n isalignment of the parachute and SRB or at a frequency of 
one-half the SRB natural frequency. 

No significant alteration of the non-elastic response characteristics of the 
SRB/ Parachute combination was seen when the elastic model was employed. 
This is not unexpected since the variations in the suspension lines and 
riser lengths are quite small compared to their steady state lengths. 

The Euler angle and angle of attack responses of the SRB/ Main Parachute 
combination for pendulum and scissors initial conditions are shown in 
Figures 37-40. 


LINEARIZATION OF THE NOMINAL DESCENT PHASE 

The linearization techniques described in Section II were applied to a variety 
of cases to obtain Root Locus Plots. Using the frozen point spectrum analy- 
technique as described in Reference 11, the eigenvalue time histories 
for both perdulum and scissors type initial conditions are shown in Figures 
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Figure 35. Main Parachute - SKB Response to Scissors 
Initial Conditions in a Non Steady A ir Mass 
(X-SRB, 2-sec Intervals) 
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Figure 36. Main Parachute * SRB Angle of Attack Response 
to Scissors Initial Conditions in a Non Steady 
Air Mass (X-SRB, 2- sec Intervals) 
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Figure 37. Main Parachute - SRB Response to a 20-deg 

Pendulum Initial Condition, Elastic Suspension 
System (X-SRB, 2-sec Intervals) 




Figure 38. Main Parachute - SRB Angle of Attack Response 
to a 20- deg Pendulum Initial Condition, Elastic 
Suspension System (X-SRB, 2-sec Intervals) 
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Figure 39. Main Parachute - SRB Response to Sci 
Initial Conditions, Elastic Suspension 
(9 1 « -20 deg, 0 3 * +20 deg Initially) 
(X-SRB, 2- sec Intervals) 






41 to 44. As expected, the eigenvalues describing the fundamental oscilla- 
tory modes cover a wider range for scissors initial conditions before 
settling to near the eigenvalue resulting from a vertical steady descent. 

The long period modes (parachute) are stable in all cases. The short per- 
iod mode describing the riser is stable with very slight damping. The SRB 
short period mode, while unstable in the initial transient response to 
large scissors initial conditions, is after a short time stable and damped. 

In viewing the eigenvalue time histories, it is important to ecall some 
important features of the linearization technique used. 

• The exact nonlinear state of the entire system is the 
reference state about which the linearisation routine works. 

• The roots to the characteristic polynomial (the eigenvalues) 
are determined from manipulation of the matrix of first 
partial derivatives which is found by applying small distur- 
bances to each of die nonlinear state variables about the 
reference state. 

• The resulting eigenvalues can each be related to a funda- 
mental oscillatory mode of one of the state variables. 

• The location of a single eigenvalue in the complex plane 
represents the local stability characteristics of the state 
variable it is associated with with respect to the exact non- 
linear condition of that state variable from which the eigen- 
value was calculated. 

• The overall stability of the entire system is a function of 
the interaction of all the nonlinear motions. 


Stability with Respect to Non-Steady Air Mass 

Figures 45 to 48 show eigenvalue time histories for the SRB/ Main configura- 
tion with no initial disturbance and a +20 deg pendulum disturbance. In a 
non steady air ma 38 the stability of the response indicated by the eigen- 
values is demonstrated through the transient response and the first gust at 
15 sec. 


Stability with Respect to Elasticity 

The eigenvalue time histories for the principal oscillatory modes of the 
SRB /Main Parachute combination with elastic suspension system when a 
pendulum initial disturbance is applied are shown in Figures 49 and 50. 
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Figure 41. Short Period Eigenvalue Time Histories for the SRB/Main 
Parachute Configuration, Pendulum Initial Conditions, 

(e. = +20 deg, 0 ., = +20 deg) 
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gure 42. Long Period Eigenvalue Time-History for the SRB / 
Main Parachute Configuration, Pendulum Initial 
Conditions (9 1 = 20 deg, 0, * 20 deg) 
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Figure 43. Short Period Eigenvalue, Time-Histories for the SRB / 
Main Parachute Confip:ration, Scissors Initial 
Conditions (0- - -20 deg, 0~ = +20 deg) 
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Figure 44. Long Period Eigenvalue Time H ;'ory for the SRB/Main 
Parachute Configuration, Sc is so. a Initial Conditions 
(0^ =-20 deg, 0j = +20 deg) 




Steady Air Mass 
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-Steady Airmass (0. = 20 deg, 0« « 20 deg) 




Parachute Configuration, Pendulum Initial Conditions 
Non- Steady Airmass (0. = +20 deg, 0« = +20 deg) 
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Parachute Combination, Pendulum Initial Conditions, 
Elastic Suspension System (0- = 0, * +20 deg Initially) 
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Figure 50. Long Period Eigenvalue Time History for the SRB/Main 
Combination, Pendulum Initial Conditions, Elastic 
Suspension System { 6 . = = +20 deg Initially) 


When compared with Figures 41 and 42, no degradation of stability because 
of elasticity is seen. 


LIMIT CYCLE RESPONSES 

Throughout the investigations of this particular recovery system, special 
attention was paid to the possible occurrence of limit cycles. In no case 
treated has a limit cycle been observed or eigenvalues calculated which 
would indicate long-term undamped oscillatory motion of any component of 
the system. 


CONCLUSIONS 


In all cases tested on the nonlinear computer simulation program, the 
recovery configurations were stable. The cases tested represent the full 
range of expected disturbances. From the 6000-ft altitude at which the 
main parachutes are deployed, the recovery system would reach a vertical 
descent attitude if it were not for the wind. The response to the wind 
causes gliding down wind. The trajectory is determined by the vertical 
descent rate and the wind speed. 

Although additional dynamics are induced by the elasticity of the suspension 
system, the overall response is not adversely affected. Large spring con- 
stants should be used to avoid sling-shot effects during transient periods of 
response. 


RECOMMENDATIONS 


The development of the present math model and computer simulation paves 
the way for useful extensions and generalizations of the analysis to provide 
a more complete and realistic representation of the entire recovery pro- 
cess including the Opening Dynamics phase. 


INCORPORATION OF PARACHUTE OPENING 
DYNAMICS IN THE MATH MODEL 

An important consideration in the overall dynamics of the parachute recov- 
ery process is the deployment and inflation of the parachute, the process 
referred to in the literature as Opening Dynamics. 

An opening dynamics analysis would establish the most realistic initial 
conditions possible by including the inflation process of the deceleration 
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system. The period in the descent phase between drogue stabilization of 
the SRB and fully inflated main parachutes sees the speed of the SRB drop 
dramatically. The dynamics of this period as described by an opening 
dynamics model would furnish more accurate initial conditions for the final 
descent and water impact. There are several Opening Dynamics theories 
which employ such concepts as dimension less parachute filling time, 
canopy volume as a function of filling time, drag areas and drag coefficient 
as functions of filling time, etc. Factors affecting the dynamics of the 
opening parachute are the canopy mass, suspension tine mass, included 
and apparent masses, and moments of inertia both real and apparent of the 
inflating canopy. Experimental data have been collected and empirical 
models nave been developed. 

It appears, therefore, very desirable to add the parachute Opening Dynamics 
to the computer simulation model based on state of the art models and in- 
cluding snatch force and opening shock calculations for the inflating para- 
chute through reefed stages to steady state. 


RELAXATION OF GEOMETRIC CONSTRAINTS 

By relaxing geometric axial symmetry constraints of the present math 
model, greater realism and additional flexibility would be obtained for use 
in stability and design analysis. 

If one allows off-axis of symmetry attach points on the SRB and the con- 
fluence point, then individual suspension line stretch and stretch rates must 
be accounted for. 

Another possible generalization would consider the parachute and/or the 
SRB to have a plane of symmetry instead of an axis of symmetry. Such a 
generalization increases the complexity of the analysis and permits the 
consideration of "gliding" decelerators and/or finned SRBs. 
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APPENDIX A 

DOCUMENTATION OF THE PARACHUTE 
DYNAMICS AND STABILITY ANALYSIS 
PROGRAMMING SYSTEM 


Computer programs describing the descent dynamics and stability analysis 
of a parachute payload system are described. 

The overall program is called CHUTER. The programs are developed in 
FORTRAN IV programming language. There are several running mode 
options. The basic running mode (no supplementary options employed) is 
simply a nonlinear dynamic simulation. Three supplementary options can 
be attached to the basic running mode. 

• Elasticity. The use of the elastic option causes the riser 
and suspension lines to become dynamically elastic and 
the nonlinear simulation to reflect the influence of the 
additional dynamics. 

• Non Steady Air Mass. The use of the non steady air mass 
option enables the subroutines describing wind and gust 
conditions to be imposed on the descending recovery sys- 
tem. The aerodynamic effects of the imposed non steady 
air mass are then accounted for. 

• A third supplementary option enables the linearization 
subroutines to be incorporated in the analysis. Their use 
causes the nonlinear equations of motion to be linearized 
at intervals in time using as a reference state the exact 
nonlinear state at the particular time. Eigenvalues for 
the linearized equations of motion are determined. 


OVERALL PROGRAM ORGANIZATION 

The overall organization finds the main program directing and controlling 
the subsequent operation of the several subroutines as well as data input 
functions. The overall organization is diagrammed in Figure A1 showing 
the subroutines and available analysis options. 

The principal variables describing the state of the system are contained in 
the "Y-array" and are passed through the various subroutines in the common 
block: 


COMMON/AAB/Y(33). 
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MAIN 


Simulation 



Figure Al. Overall 
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The time rates of change of the state variables are contained in the "D>array" 
element having the same index and are passed through the various subroutines 
in the common block: 


COMMON/AAC/D(30) 

The principal variables are listed in Table Al. Nearly ail other parameters 
and variables and constants which are required by more than one subroutine 
are passed through a series of common blocks containing related arguments. 


CHUTER INPUT/OUTPUT 
Input Description 

An input card deck of 14 cards provides the required information for initial- 
ization and control. The input data deck is described in Table A2. 


Output Description 

There are two forms of information output from CHUTER. The line printer 
output provides detailed information on the nonlinear simulation at chosen 
time points along the trajectory, the interval being DTP. When the lineari- 
zation option is employed, the eigenvalues of the linearized system a~e 
printed for the points along the trajectory at which the nonlinear m is 
linearized. 

A plotting subroutine is included which charts information generated by the 
nonlinear simulation subroutines. Additional charts are drawn if the elastic 
or non steady air mass option is employed. 

The line printer output during nonlinear simulation consists of groups of 
four lines each corresponding to the time printed at the left of the page. 

Each page is headed by column labels. 

When the linearization routines are employed, the eigenvalues at the se- 
lected linearization points are stored until a single page can be printed with 
the eigenvalues for the previous five linearized points. 

For each new run a run title page is printed listing the supplementary op- 
tions employed, and a data deck reproduction is made for reference. An 
illustration is drawn on which the principal system initial geometric para- 
meters are noted. 

The line printer output continues through the maximum simulation time or 
water impact. The exact state at that point is printed. 
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TABLE A1 - DEFINITION OF PRINCIPAL VARIABLES 
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TABLE A2 - TYPICAL DATA CARD INPUT DECK 



3 





TABLE A1 - TYPICAL DATA CARD INPUT DECK 
(CONTINUED) 


Data Description (All Data is Floating Point) 

Card 1. Format (2F 8. 0) 

Variable 

Units 

Definition 

Y (30) 

ft 

Initial altitude 

HDOT 

fps 

Rate of Descent 

Card 2 Format (6F 8. 0) 

Variable 

Units 

Definition 

D3 

ft 

SRB Diameter 

L3 

ft 

SRB CM Location 

L3T 

ft 

SRB Total length 

L4 

ft 

SRB CP Location 

M3 

Slugs 

SRB Mass 

S3 

ft2 

SRB Cross Section Area 

Card 3 Format (3E 10. 3) 

Variable 

Units 

Definition 

DOC 3 

slug ft 2 

SRB Inertia about its X axis 

IYY 3 

slug ft2 

SRB Inertia about its Y axis 

IZZ 3 

slug ft 2 

SRB Inertia about its Z axis 

Card 4 Format (8F 8. 0) 

Variable 

Units 

Definition 

BCN (Array) 


Constants in the polynomial 
describing the normal force 
coefficient of the SRB 

Card 5 Format (9F 8. 0) 

Variable 

Units 

Definition 

BCT (Array) 

m mm 

Constants in the polynomial 
describing the tangent force 
coefficient of the SRB 
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TABLE A1 - TYPICAL DATA CARD INPUT DECK 
(CONTINUED) 

Data Description 

Card 6 Format (9F 8. 0) 

Variable Units Definition 

BCM (Array) — Constants in the polynomial 

describing the moment coeffi- 
cient of the SRB 


Card 7 Format (9F 8. 0) 

Variable 

DO 
LI 

LSO 

M 
MC 
ML 
LCM 

SI 

CLUST 


Card 8 Format (9F 8. 0) 

Variable 
ACN (Array) 


ACT (Array) 


Card 9 Format (9F 8. 0) 

Variable Units Description 

ACM (Array) — Constants in the polynomial 

describing the parachute 
moment coefficient 


Description 


Constants in the polynomial 
describing the parachute normal 
force coefficient 
Constants in the polynomial 
describing the parachute tangent 
force coefficient 


Units 


Description 


Parachute nominal diameter 
Length from confluence point to 
parachute CP 

Initial value of suspension line 
length 

Number of suspension lines 
Mass of canopy 
Mass of lines 

Initial guess at parachute CM 
distance from confluence point 
Nominal parachute area 
Number of chutes in cluster 




TABLE A1 - TYPICAL DATA CARD INPUT DECK 
(CONCLUDED) 


Data Description 


Card 10 

(Format (IF 8.0) 



Variable 


Units 

Description 

L20 


ft 

Nominal riser length 

Card 11 

Format (6F 8. 0) 



Variables 


Units 

Description 

Y (4) 


deg/ sec 

Initial PI 

Y (5) 


deg/ sec 

Initial Q1 

Y (6) 


deg/ sec 

Initial R1 

Y (7) 


deg 

Initial el 

Y (8) 


deg 

Initial 01 

Y (9) 


deg 

Initial 41 

Card 12 

Format (6F 8. 0) 



Variables 


Units 

Description 

Y (13) 


deg/sec 

Initial P3 

Y (14) 


deg/ sec 

Initial Q3 

Y (15) 


deg/sec 

Initial R3 

Y (16) 


deg 

Initial <j3 

Y (17) 


deg 

Initial 03 

Y (18) 


deg 

Initial 43 

Card 13 

Format (IF 8. 0) 



Variable 


Unit 

Description 

TMAX 


Sec 

Maximum time 

Card 14 

Format (3F 8. 0) 



Variable 


Units 

Description 

YWIND 


— 

Wind option, 1-yes, 0-no 

YELAST 


— 

Elastic option, 1-yes, 0-no 

/LIN 


m m 

Linearization option, l-yes 
0-no 





A typical page showing information on the nonlinear simulation is shown in 
Figure A2, and a page showing the eigenvalues at selected points is shown 
as Figure A3. 


PROGRAM DESCRIPTION 
Main Program 

CHUTER is a series of subroutines whose operation is controlled by the 
MAIN program to provide nonlinear and linear analysis. The MAIN pro- 
gram is diagrammed in Figure A4 and a source listing is presented in 
Figure A5. 

The MAIN program is broken down into three parts. The first is input and 
establishes constants and control variables. The second segment initial- 
izes the elastic variables, sets angles and angular rates to units of radians, 
sets the initial velocities in the body fixed coordinates, and establishes the 
direction cosines matrix corresponding to the initial conditions. Finally, 
the third segment is a high-frequency loop which runs the nonlinear simu- 
lation. 

The high-frequency loop is initially entered with mode and time = 0, which 
causes the initial conditions to be output by subroutine PRINT. Successive 
passes through the loop increase the MODE to its nominal value of 4 or 5 
depending on whether subroutine PRECOR is about to predict or about to 
correct. 

Elasticity initial conditions (i. e. , riser and parachute center of mass 
lengths) are updated through time - 0. 25, at which point the numerical 
determination of elastic rates begins. 

Time = 0. 25 is an arbitrary but convenient time greater than time = 0 since 
at time = 0 the elastic elements are unstressed. 

There are four normal exits from the high-frequency loop. After the print 
time interval DTP the loop is exited by a call to subroutine PRINT. The 
•ccond normal exit occurs when a water impact occurs. This is sensed by 
comparing the altitude with length from the SRB center of mass to the 
engine end. The third normal exit occurs when the simulation time exceeds 
TMAX. The fourth normal exit occurs when a point in time is reached 
about which a linearized solution is to be found. 

The subroutines used with CHUTER are listed in Table A3. 
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Figure A2. Typical Nonlinear Simulation Output Page 
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Figure A3. Typical Eigenvalue Output Page 








Figure A4. CHUTEB Main Program Flow Diagram (Continued) 













Figure A4. CHUTER Main Program Flow Diagram (Concluded) 
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PROGRAM CHUT ER C ! NPUT » OUTPUT » TAPES* l NPUT * T APE6-0UTPUT .TAPE?) 

r 

r OEFINITIONS op PRINCIPAL VARIARLFe 

r vm • ui y c i o > = it? 

r YC ? ) » VI YC?0) « V2 

C YC A ) « W1 Y C ?1 ) « V? 

r yca) * pi yc??) * p? 

c YCA) « 01 Y ( 11 ) * O? 

C Y (A) a 0) Vtm « P? 

YC7) a PMlI Y C ?S ) « PHt? 

Y( A ) * THPTal YC?M « THETA? 

YC ° ) * PSll YC?7) « P?I2 

YUO) a U3 YC?M * SRft CG POSITION. EARTH X AXIS 

Y( 1 1 ) a V? YC?R) a SRR CG POSITION. EARTH Y AXIS 

YU?) a WA YC AO) a SR» CG POSITION. EARTH ? AXIS 

Yl 1 1 ) a PA YCA1) = PARACHUTE CG POSITION. EARTH X AXIS 

YU4) * OA YC A?) = PARACHUTF CG POSITION. EARTH Y AXIS 

Yl 15) a RA YIAA) a PARACHUTE CG POSITION. E AR TH Z AXIS 

YC1M * PMI3 
YC1 7) a THETAS 
YC1«) a DSIS 


1 PARACHUTE 

- ? O t.SED 

A PAYLCiACi 


COMPLETE LIST OF all COMMON ELOCXS 

COMMON /A A A /ACT C9 ) » ACNC 9) . ACMC R) «RCT C « ) .RCN I 0 ) »RCMC 0 ) 

COMMON/A A°/ YC3M 
rOMMOM A A AC / n C “*0 ) 

COMMOM/AAf'/S C *«B.M «f>SCf ) »TC A.M 
'■OMMON/aae/ AAC9) »AC|Q) ,ACC9) 

C0MM0N/AAF/CN1.CT1 .CMS »CTA * ALPHA1 .ALPHAS . APT A 1 .EfTAA .GAMMA 

COMMOM/AAEF/CMl.CMA.ALCM 

COMMON/AAG/L ? »L?O0T .L2000T .LCM.LCMOOT .LCMDOT 
C0MM0N/AAH/C1 .C3.E2.L3.RA0.L1.LA.CF1.CF3.S1 .S3 

C omvov/AAHH/ii 

COMVON/AAJ/VOnp 

C0MM0M/AAX/AL1 »AL?.AL3.ALA 

COMMON/ A AL/F 1X.F1Y.F1Z. FAX. FAY. FS/ 

C0MM0N/AAM/M1X»M1Y.miZ.MAX.M3Y.MAZ 

C0MMON/AAN/A5.A7, AA.A9. A 10. All *A1 2.A1A.A14. A1 S, A 1ft, A 17. AIR. A1R.A20 
CO MM0A’/AANN/A?1 

COmmom/AAO/OI "TJME.YPLAST.ETIMP 
rOMMOH/AAP/YwiNO»VV*IA'O.VG"eT. w IGi' 

COMMON/ A AO/ I XXI .IYY1,IZZ1»I*X3.IYY3.I?Z? 
COMMON/AAQO/IXZ1.1YX1.IZY1.IXZ3.IYX3.IZY3 

COMMON /AAR/KLS. KR.MC. ML .Ml .MIA .MP.OO.RO.RHO.L20.LCMO.LSO.M.M3 

common/aarr/ls.lcp 

COMMON/A AT/T IMF 

COMMON/ A AW/NO°S 

COMMON/PLTR/XXI40?) .THF1C402) »THF3C40?).AP1C4G2) .AP3C402* .ALTC40?) 


Figure A5. Main Program Source Listing 
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1 *RNO(402) »EOR(40?) »RL(40?) *CL(40>) ,Wfi(40?) 

C0MM0N/X0R0S/SUMMA1 *SUMMA?»TOTAL»AVERAl »AVFRA2»OVERAl *OVERA2»DT 
RFAL TXXAl*tXXA3»TXXl»IXxio*IXX3»!XLl»IXZ3 
RFAL !YYA1*!YYA3*IYY1*IYY10*IYY3,IYX1,IYX3 
RFAL IZZA1,TZZA3.!ZZ1.TZZ10,IZZ3*IZY1,IZY3 

REAL KLS*KR*LCM*LCMDOT ,LCMDDT »LCMO,LCP*L5»LSO*Ll *L2»L2DOT *L2DOOT 
REAL LZ0»L3*L3T»L4»MC»ML*M1»M1A»M1A»M1Y,M1Z*M3X»M3Y»M3Z*NP*M3*M 
DT*1,0 

dtx*7,o*pt 
1 f,oto ’ 

3 '••RITr (6*430) 

Of) 7 l*NOPS»400 
THFI(J) ■ THF1 (NO«S) 

TMF3(T) m THF* (NO°S ) 

ADI (1) * API (MO«S) 

ADI (1) * ADO (NODS) 

PMF- (!) * R*T, (NODS) 

ALT (!) « ALT (NODS) 

FDD (!) * FOR (NODS) 

RL (!) « RL (NO«S) 

CL (!) * CL (MO«S) 

7 WO (!) » WO (NO«.S) 

FALL PIFTIIR ( YWINO,YDLAST) 

3 CONTINUE 

sijMMAl* SI JMMA 7* TOTAL* o.o 
r 7FR0 Y AMO D ARRAVS 

or* ID ! a 1,30 

o( T ) » 0.0 

10 Y(!) * 0,0 

C RE AO initial ALTITUDE AND RATE OF DESCENT — RFCALL THAT THE +2 AXIS OF 

c thf earth is directfo downward so that the initial altitude is negativf 

DPao (3,410) Y ( 30 ) *HOOT 
r DT A O «n PHYSICAL DIM^N'TONS 

pf AO (3*410) 03*L3,L3T»LA,M3»S3 

F DFAD ten INCRTIAL CHAR ArTFQ I ST ICS 
RF AD (3,420) IXX3»IYYX*IZZ3 
F RFAO AFRoDYNAHIF FHARAFTFR! ST *CS 

RFAO (3*410) (40M(I),I«1»A> 

dfao (*,410) (PFT(!)*!*1»D) 

dpap( 3,410) (BfM( !)«!*) ,0 ) 
r r>FAO DAPArwijTF PHYSICAL dimfnsIONS 
C CLUST IS THF NUMBER OF FHUTFS IN THF CLUSTER 

RFAO (3*410) D0,L1,LS0,W,mC.ML,LCW0*S1»CLUST 
C READ PARACHUTE AERODYNAMIC CHARACTERISTIC EQUATIONS C0FFFIC1FNTS 
RF AD (3*410) ( ACN( I )*I*l»3)*tACT(l)»I*l»6) 

DFAO (3*410) (AFM(U , !«),D) 
r Rfa" SYSTEM CHARArTrRlSTIFs 
RFAO (3.410) LZO 
C READ INITIAL CONDITIONS 

RFAO (3*410) (Y(t)» 1*4*9) 

DFAD (3,410) !Y( T ) *1*1 3,1 A) 

r RFAO TMF TIME AT WHICH THF ANALYSIS MUST STOP 
DPAO (3*410) TMAX 

C READ WINO, FLASTIFITY, AND LINEARIZATION OPTION CONTROLS — 0 FOR NO AND 1 FOR 


Figure A5. Main Program Source Listing (Continued) 
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C VFS 

READ (5*410) YWINO.YELAST.YLIN 
C DATA CARD LIST 

WRITF (6*550) 

WRITF (6,j55) 

IF (YHIND.NF. 0*0) WRITF (6*556) 

IF(VFLAST.NF.0.0)V*RITF (6*557) 

!F(YLTN *MF.0«0)WRITF (6*65R) 

IF| YLlN*FO.O.n,aND.YFLAST.FO,^ # 0«AND«YW!ND,Ff),r>,o) yRJTF (6*550) 

WRITF (6*560) * ' 

WRITF (6*430) Y(30),HD0T 

WRITF (6*430 ) D3»L3*> 3T*L4*M3*S3 

WRITF (6*440) IXX3*IYY3*IZZ3 

WRITF (6*430) (BCN( 1 ) * 1*1 »ft) 

WRITF (6*430) (RCT( I )*T«1*0) 

WOITF(6*430) C*CM(I),I«i*9) 

WR I TF (6*430) no»Ll»LSO,'SMC,VL*LCMO*Sl*CUiST 
WRITF (6:430 ) ( ACM I ) * T«) *3 ) . ( ACT( I ) » 1*1 *6) 

WRITF (6*430) ( ACM) I ) ,!«1»9) 

WRITF (6*430) L70 

WRITF (6,4?;)) ( Y ( I ) * I *4 *0 ) 

WRITF (6*430) (Y(T)*ts)3*]ft) 

WRITF (6*430) TWAX 
WRITF (6*430) YwlND*YEL AST »YLIN 
C CONSTANTS 

DTO . 0*75 
ftjvf « 0.75 
O » V,17 

J • 0 

LCM » LCWO 
LFP « 0*163*00 
LS • LSO 
l 7 . » 170 

moof ■ o 
mo « Mr *r lust 

WL « ML*CLUST 
PI ■ MC*WL 
N • 30 

OLOTIMF ■ 0,?5 
PAO ■ 57* ?057P 
R C « 1,0 

RO ■ 0*36*00 

SI ■ S1*CLUST 

TIMF « 0*0 
VWt mo « VO* 1ST ■ 0,0 
WTOli * 0,0 

r aliow p t sfr with thf RARACHurr axis of symmftry 

OO 15 T* 77, Z7 
15 Y( I ) > Y(!-18) 

f PLASTICITY INITIAL CONDITIONS 
LCWOOT * 0,0 
LCMODT « 0.0 
L’OOT • 0,0 
LTOOOT 0,0 


Figure A5. Main Program Source Listing (Continued) 
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AL1 ■ L> 

* L200T 
AL3 ■ L r»* 

AL4 « lcnooT 

C *L*STIC ro'FFtrp'MTS 

if (''*last*rt.o.o) goto I 7 

<R * l.o*?n 

<L5 • 1.0*70 

ROTO |& 

13 CONTIAUr 

<R « V3»R*)S./L30 

icls • w3*r*15./*/lso 

14 rORTTVUF 
CALL TOR AO 
CALL SR**t»> 

fall i«iv*l« funou 

ROT* 31 

’0 calc oirrps 

31 fORTIRUF 
C UPOATF L7.LC* 

IFITTNF.RT.O.IS) o-OTO a? 

31 M.1 * L7 

AVI * LC 1 ' 

37 CONTIN'IF 
C WIRO CALLS 

IF (VWINO, *0.0.0) ROTO 17 
CALL Wll»o 

IF (TIM* .*0. 0.0) r.OTO 1*. 

C RUST CALLS — 0*9100 » 15 **C 

IFfTIMF-WTIMF-15.0) 17.16.16 

16 CALL OUST 

WTIRF « TIMC 

17 •'ONTIMUF 

WIRII * VWINO ♦ WRUST 
CALL CHUT* (CLUST) 

CALL COFFTS IRO 
call FORCFS 
CALL MOMFNTS (RHO) 

CALL OIFFOR 

IF (TT*F,nc.o, 0) r-OTO no 
C ILLUSTRATION 

WRITF 16.1000 JCLUST. 01. SI. MC.LCP.R0.ML.LCM.LS.L2.03. S3. LIT. L4.L3* 
1*3 

GOTO 120 

R9 CALL PRFCOR (N.^.M) 

C U2.V7.W7 

HO 1 75 I « 1*3 

125 YIlR*t>« (Yd )*Y<S>*LCM>*(«M1.1 ♦ !)* PI2.1.1 )*P( 1 .?♦ !)•«( 7.7.1 1 
1 +6(1 .3. !)**»(?. 3.1 >)*(Y(7)-YU)«LCM)»(B(1.1 .I)*B(?.?.l) 

1 ♦B(1.2.1)*P(2.2.?)*R(1.3.I )#P(?,3,?))-.Y(3)*(B(1»1 »t 1 

1 *R(2.1.3)*P(l»7.I)*P(?,7»3>*n{l»3.I)*R(?,7,3)) 

IF (Y(30) ,RF. -75.0) ROTO 170 
110 IF CTIMF-TPRINT-0TP170 .170.170 
120 CALL PRINT 


Figure A 5. Main Program Source Listing (Continued) 



T»RTMT ■ TfNF 
IF (YLIN.F 0 * 0 . 0 ) r,nrn no 

POL I S*I TR)t<*TtMP 

IF( JPOLIF.FO.TIMC j.A>m.(TIvE.NF.P.O> )CALL OPR I VF I POL I S*RC »RHO *rLl|F 
IT) 

130 CO NTTWUF 

!FIvn0),6F,-?5.0) ROTO 5 
IF IT IMF »6F. TMAXI GOTO 5 
IF ITfMFl gO,R9,in 
410 FORMAT IIOFU.O) 

420 FORMA TO?) 0,3 ) 

430 FORMAT I10F9.3) 

4*0 FORMAT !3F12#3I 

450 FORMAT ! ?0X* 3 1HWATFR IMPACT OR MAX T I MF — cnoj 
550 FORMAT ( 1 Hi ) 

555 FORMAT! 10X*41hPARACM!)TF DYNAMICS A NO STABILITY 4NALVSF S* //, 1 ?x . 
125MNOM L INFAR SIMULATION wjTM,/) 

556 FORMAT I 1 ?X » 1 9HN0N STFAOY AIR MASS* / ) 

55T FORMAT »1 2X »34HFLASTIC RISCR AND SUSPFNSION LINFS,/) 

558 FORMAT II 2X .42HNUMFR ICAL L I NFAR 1Z AT I ON AT «FLFCTFD POINTS,/) 

S5R FORMAT! nx • I OMNO OPTION?,//) 

560 FORMAT t’X* 1PHDATA CAROS AS RFAO,//) 

1000 FORMAT! 1H1»37X»F?*0»1 4H CH'.lTF CUISTFR ,7X,?HCCCCrCC* 11X »4HP1 «,C7,1 
I*3H FT«/,57X.1HC»13X*1HC»TX»4HS1 =*F?,1,6H s C F T»/,?5X * 1HC« 1 7X , 
?1HC »/*54X*lHC »19X*1HC »4X »4HMC * »F7, 1 *6H SLUGS, /»53X *1HC*21X* lHC»/» 
352X*1HC*7X*SHCP 0.UX.7HC LCP =*F7.1,3H FT */*S1X,1hC.25X*1mc,/* 
451X,lHC,?5X»lMC./,57X,?5MrcCCCC0CCCCrrCCCCrrCCCCCC,/,5ZX,lHL,7'»x,7 

5HL RO *»F7» 1 ,3M FT »//»53X » l^L* 7 1 X * IHL */» 7ox ,&HML <xnGS, 

6/*54X*lHL»19X»lHL»/,ftOX,5HC |, 1 0*’ 4X *4MLC*= *F7, 1 ,->H cy , /,ssx , 1 u L , 

717X,1HL»//,56X»1HL»15X»1HL»//*57X»1HL»1?X» , ML Ls=,F 7,1,-»H FT,//, 
858X»lHL*nX*lHL»//»59X,lHL*9X,lHL*//*60X.lHL*7X,!HL*//*61X»!HL, 

95X, 1HL»//*62X »1HL» , X»1HL*//*63X*1HL L »//»3PX » 5HR I ?FR * ?1 X »fiHP L? 

1“*F7.1*?H FT ♦/*«! 64X»IfR»/) #36X»? 'HSOL IP ROCXft 3 00ST C, ?*6X* l >- , f , 

14H03 = »F? , 1 « 3m FT,/*63x»Tpo 0,8X,4H6? =,F7.1,fH 'Q Fl,/, 6 ->X, 

25H« B»/,62X.5M« ° »7X *AHL3T« ,F7, 1 , *H FT , / »R ( 62X * sun "./), 

356X ,9HCP B 0 B»7X»4ML4 »,F7,1,3m FI*/»67X,5HF P,/,5AX, 

49HCM B 0 B*7X*4HL3 *.f7.1,3M FT,/,6?X*5M? p,/,6?X,5hR P,7x, 

54HM3 «»F7.1»6H SLUGS*/ ,5 ( 62X ,*H p p. t »62 X»fhprpor,/) 

pro 


Figure A5. Main Program Source Listing (Concluded) 
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TABLE A3 - BASIC SUBROUTINES USED IN PROGRAM CHUTER 


Symbols 
Table No. 

o h n m t> to - . 

tf) <0 fCDftHrtHHrtH • ! 

< < < <<<<<<<<<• ' 

• 

o 

a * 
• 3 

2*? 

fa 

< < < <<<<<<<<<<< 

Diagram 
Figure No. 

Afl 

A8 

A10 

A12 
A14 
A 10 
A18 
A20 
A22 
A24 
A26 
A28 
A30 
A32 

r" 

Description 

Differential Equations of Motion 

Parachute Geometry, Inertia 

Parachute and Payload Aerodynamic Coefficients 

Aerodynamic Forces 

External Moments 

Direction Cosines 

Direction Cosines Rates 

Predictor- Corrector Integrator 

Mean Wind Profile 

Gust Envelope 

Elastic Rates 

Output 

Constants 

Initialize Velocities 

Subroutine 

DIFEQN 

CHUTE 

COEFTS 

FORCES 

MOMENTS 

DIRCOS 

DBDT 

PRECOR 

WIND 

GUST 

ELASTIC 

PRINT 

CONST 

INVELO 
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Basic Subroutines 


The basic subroutines are those which describe the aerodynamics, the dy- 
namics, or the kinematics of the nonlinear simulation, the nonsteady i c 
mass models, and techniques used in the linearization of the equations jf 
motion. All the other subroutines are manipulatory in nature and hence are 
termed auxiliary subroutines. 

Subroutine DIFEQN -- Subroutine DIFEQN implements the system of differ- 
ential equations L Equations (16) to (39), (52) to (69)J. The time derivatives 
of each of the state variables and the riser force are calculated. Moments 
about the body fixed axes for the parachute and SRB are updated due to tne 
change in riser force. During the Runge-Kutta initialization steps anu the 
predictor step of subroutine PRECOR, the section of DIFEQN containing the 
equations coupled by the riser constraint is looped through four times to 
ensure that the influence of the coupled terms is uniform. Subroutine 
DIFEQN is diagramed in Figure A6 and a source listing is presented in 
Figure A7. Table A4 presents a list of symbols for DIFEQN. 

Subroutine CHUTE -- Subroutine CHUTE computes the geometric and initial 
characteristics of the parachute as a function of time. Also calculated is 
the air density as a function of altitude. 

The parameter CLUST, passed in calls to CHUTE, represents the number 
of chutes in the cluster. As all the input data were for a single chute, the 
mass and inertia are multiplied by CLUST to form the mass and inertial 
characteristics of the single chute equivalence to the cluster. 

The parachute center of mass location is calculated as a function of the 
canopy mass, the suspension line mass, and the mass of the air included 
in the canopy. 

Finally, when the elasticity option is employed, the ELASTIC subroutine is 
called to compute the rates of change of the lengths of the elastic elements. 

Subroutine CHUTE is diagramed in Figure A8 and principal variables are 
defined in Table A5. A listing of CHUTE is given in Figure A9. 

Subroutine COEFTS -- Aerodynamic coefficients are calculated for the 
normal and tangential forces and the moments ou the parachute and payload. 

The coefficients for normal force and moments are calculated as a function 
of the angle of attack, a, using the polynomial form. 

C| or + Cg cr^ + Cj + + Cg + Cg 
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Figure A6. Subroutine DIFEQN Flow Diagram 













SUP ROOT INF niFEON 

C SUBROUT INF DIFFON CONTAINS THE STSTFN OF NONLINEAR DIFFERENTIAL EQUATIONS C 

C IT RETURNS THF SLOPES OF THF FUNCTIONS FOR THF TIMF AT WHICH IT IS CALLEO C 

COMMON/ AAP/Y 03) 

C0MMON/AAC/D(30) 

COMMON/AAO/B ( 3*3*3 ) *BS(6 ) »T( 3*6 ) 

COMMON/AAF/ AA(9)*Aft(9)*AC(9) 

COMMON/AAG/L 2 *L 200T *L2DOOT »L<M*tCM00T *LC MOOT 

COMMON/AAH/CltC3*F?»L3»RAO»LltLA»CFltCF3*SltS3 

COMMON/AAJ/MOOF 

COMMON/AAL/FlX»FlYtF12»F3X»F3Y*F3Z 
COMMON/A AM/M 1 X »M1 Y tMl 2 »M3X »M3 Y »M3Z 

C0MM0N/AAN/A6*A7,A8.A9»Al0*All*A12»A13.AlA,AlS*A16,Al7,Alft*Al9,A?0 

C0MM0N/AANN/A21 

C0MM0N/AAQ/rxXl.!YY).t2Zl»!XX3.IVY3»IZ23 
COMMON/ AAOO/ IX21»IVX1»I2Y1»IX23»IYX3*I2Y3 

COMMON/AAR /KLS»KR»MC*ML»M1.M1A*MP, DO. R0.RM0»L?0*LCM0*LS0»M,M3 
COMMON/ A ARR/LS*LCP 

REAL I XX| • I YV1 . 1 2Z1 * I XX3 •! YV3 • 1 223 

REAL IXZltIYXl.|2Yl»IXZ3»tYX3»I2Y3 

REAL Ml#M3»NP,L3»XLS*KR»MC*ML»MlA.L2DOTtL20DOT»LCM,LCMOOT»LCMOf>T 
REAL L1*L3 T»l4*L?0*LCM0*LS0 
RFAL 12*LCP*LS 

REAL Ml X *M1 Y *M|2 *M3X *M3 Y »M3Z *M 
C CONSTANTS IN THE OIFFFRFNTIAL FQUATIONS 


G 

* 

32*17 

A1 

9 

Y(10)-Y(1A)*L3 

A? 

9 

VIU )*Y(13)*L3 

Al 

9 

Y(1)*Y(5)*LCM 

A4 

9 

V(2)-Y(A)*LCM 

A* 

9 

Y ( 3 ) ALCMOOT 


CALL CONST 

CALL ONDT 

MIXF»1 

IF(MOOF*LF*A) MIKF«* 

DO BO RK*1*MIXE 
C 0W1/DT 

0(3) « ♦(l./B(l*3.3))*((AA(3)»Y(23)+Y(22>*ARUi>)*L2AL200T*2.n 

1 »AC O )-L?*(BI ?* 1»M*D( 23 ??>)-*(*♦?, 3 )»L?o 

2 POT+AA(6)*A1*A6-*AR(6)*A2+A7.*AC(6>*Y(1?)+A10-AA(9)«A3-A 

3 8-AR|9)*A4-A9-AC(9)»A5)-LCMD0T 
C RISFR FORCE 

F2 • (MP»(D(3)«-Y(2l*V(4)-Ym#Y<5)>-FlZ-m*ft(l»3»3)*G)/RS(3) 

CALL MOMFNTS (RHO) 
r 1*1 DOT 

D(l) • (FI X+M1*N( 1*1 *3 )*G*F?*BS( 1 ) ) /MP«Y (4)**Y(5)+Y(2)*Y(M 
C V1D0T 

0(2) » (F1Y-*M1*R(1»2*9)*G*F?AN$(2) )/MP«Y(l)»V(6)AV(3)*Y(4) 

C PI DOT 

0(A) • (M1X-!ZY1«V(5)*V(6))/IXX1 

C 0100T 

0(5) • (M1Y-IXZ1*Y(6)*Y(A))/IYY1 

C R1D0T 

0(6) « (M1Z-IYX1»V(A)*Y(5))/I2Z1 


Figure A7. Subroutine DIFEQN Source Listing 
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C U3D0T 

DUO) « tF3X*M3*Pt3.1 »3)*G-F2#PSt*) » /M3-Y (12|»Y(1* )4-Vtl 1 )*Y< 1* ) 
C V3DOT 

Dtll) « tF3Y4.M3#Rt7»?.3)#r,-F?*PSt*> ) /V3-YU 0)*Yt l* )+Y 1 1 ? )*Y 1 1 3 ) 
C W3DOT 

DU?) - tF3Z+M3*Pf*.3,3)*Ci-F?*PSt6> ) /M3-Y f \\ ) *Y(13 ) 4Y (10 )*Y Cl A ) 
C P3D0T 

Dt 13) * (M3X-IZY3*Y(14)*YU5))/|XX3 

C 0300T 

DU*) = IM3Y— IX?3*Y| 1 A ) *Y ( 13 ) ) / I YY3 
C R3D0T 

DUS) * tM3Z-IYX3*Ytl3)*YtlA) )/lZZ3 
CALL CONST 
C DW1/DT 

Dt 3 ) = +t1 ./Pfl.3.3) )*( t4At»)*Vt?3)+Yt??)*APU))*L?+L?DOT#9,o 

1 *ACt31-L?*tPt?»l*3)*0t?-*t-0t?.?.3)*Ot?7))-nf?,-».3)*L?f' 

2 DOT+AAt 6)*M+»A+A«'tA)**->*A7«.»f tft)*Yt 1?)+M A A t O ) * A3 —a 

3 8-A0f9)«A*-A9-art < »t*A , »)-LCM'V'T 
C DP2/DT 

Dt 22 ) = -t 1 ./Bt 2»?»?) )*tYt?3)*A*t?)+Yt?2)*A°t ?)-Dt ??) # n t ?»1 .?) 

1 +tl./L2)*t L7D0T*?. ~<*AC t ? ) - !s t ? .3 * » ) *L?DprT+AA I « ) *At *A l 

2 l+ARt5)*A?+M34.Art5)*Ytl?) + A)3-AAt A)*A-»-A)'.-ARt R>eA*-A 

3 15-Af tfl)*A3-A>1*ntl.3*?) ) ) 

C DQ2/DT 

Dt ?3 ) « + tt./Pt?.1.1))*tYt23)*AAtl>+Yt??)*A'»tl)+'>t??)* c *t?»?.1) 

1 +tl ./L?)*tL?DOT*t?.D*t/.C») ) ) )-Pt 2.3,1 )*L?PD0T4 .aa(4)*A) 

? +A16+AP t A ) * A’+A17+AC t !• i *Y 1 1 ? )4A 18-A A t T )**■>— A) D-4P 1 7) a 

3 A*-A?P-ACt ?)*AS-A?l*Pt 1.7,1))) 

80 CONTtNUF 
C PHI1D0T 

Dt 7 ) = Yt*) + 1 Y t 5 ) *SINt Yt 7) ) +YI6) *C0StYt7)) ) *TAN t Y t P ) ) 

C PHI3D0T 

Dt 16) = VU3)-MYUA)*S1NIVUA) )+VUS)*C0StYM6)) )*TANtYU7) ) 

C PWI2D0T 

Dt ?5 ) = Yt??) + tYt?7)*StN(Y(7 , >) )-f-Y ( 2* ) *COS ( Y t ?5 ) > ) »T \N t Y ( 9ft) ) 

C THFTA1D0T 

Dt 8 ) « Ytf.) *Cv ;Yt7)) -Y(6) *SIN<Yt7)) 

C THETA3D0T 

Dtl 7) « Yt 1 * ) *COS t Y 1 1 A ) )— Y 1 1 S ) *s INt Y f 1ft ) ) 

C TMFTA7D0T 

Dt 28 ) * Yt?3)*C0StYt?S) )-Yt?*)*e!Mt YUS) ) 

C PS 11 DOT 

Dt 9 ) a t Y t S ) *S!Nt Y 1 7 ) ) +Ytft) *rOSfY(7))) »SFCtYtB)) 

C PS13D0T 

DUB) = t Yt 1* )*STN t Y U6 ) )+Y(lB )*COSt Y 1 1 A ) ) )*s F CtYtl7) ) 

C PSI2D0T 

Dt 27) « tYt ?-> )#SINtYt 7? ) )+Y(?*)»rOSt YUS ) ) )*' l i<'tYt >6) ) 

C U2D0T IS Dt 19) , V2D0T IS Dt 2D) , AND W2D0T IS Dt?l) ALL UNU'-FD 
D( 19) « 0.0 

Dt 20) « 0.0 

D { ? 1 ) « 0.0 

C R?OOT 

Dt?*) « 0.0 

C X 3 FOOT 

Dt 28) • +tYU0)*Pn.l.l)+Ytl1)*nf3*7.])+YU?)ft<M3'3*l)) 

C Y3FD0T 

Dt 29) » +tY(10)«P(3»l,?)+Ynn*ot3.?.?)+Y(l?)*P{3,3,?)) 

C Z3ED0T 

Dt 30) « 4(YU0)*Pt3.1.3)4YMl)«P(3.?.?)4Yfl?)«Pt?»3»3)) 

RFTURN 

FND 


Figure A7. Subroutine DIFEQN Source Listing (Concluded) 
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TABLE A4 -LIST OF SYMBOLS FOR SUBROUTINE DIFEQN 


Quantity 

Mnemonic 

Units 

Description 


D(l> 

ft /sec 2 

Parachute CM linear accelerations 
in XYZ body fixed axes directions 

* 

V 1 

I 

ft/ sec** 

* 

*1 

1 

ft/ sec 2 


*1 

■ 

o 

rad/sec 

Parachute angular accelerations 
around XYZ body fixed axes 

«l 

D(5) 

rad/ sec 2 



D{6) 

rad/ sec 2 


*» 

D(7) 

rad/sec 

Parachute reference frame Euler 
angular rates 

h 

D(8) 

rad/ sec 


\ 

D(9) 

rad/ sec 



D(10) 

ft/ sec 2 

SRB CM linear acceleration in XYZ 
body f xed axes 

m 

V 3 

D(ll) 

ft/ sec 2 


*3 

D<12) 

ft/sec 2 



D{13) 

rad/ sec 2 

SRB angular accelerations around 
XYZ body fixed axes 

% 

D(14) 

rad/ sec 2 


*3 

D(15) 

rad/sec 2 


»3 

D(16) 

rad/ sec 

SRB reference frame Euler angle rates 

»3 

D(17) 

rad/sec 


*3 

D(18) 

rad/sec 


h 

D(22) 

rad/sec 2 

Riser angular accelerations about XY 
body fixed axes 

*2 

D(23) 

rad/ sec 2 
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TABLE A4 - LIST OF SYMBOLS FOR SUBROUTINE DIFEQN (CONCLUDED) 


Quantity 

Mnemonic 

Units 

Description 

*2 

D(25) 

rad/ sec 

Riser reference frame Euler angle 
rates 

*2 

D(26) 

rad/ sec 


• 

Y 

2 

D(27) 

rad/ sec 


X v 

E 3 

D(28) 

ft /sec 

Down range, cross range, and altitude 
rates of change of the SRB center of 
mass 

*e 3 

D(29) 

ft /sec 


Ze 3 

D(30) 

ft /sec 



F2 

lbs 

Riser force 
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Figure A8. Subroutine CHUTE Flow Diagram 
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Figure A9. Subroutine CHUTE Source Listing 



TABLE A5 - LIST OF SYMBOLS FOR SUBROUTINE CHUTE 




Units 

Description 

ALCM 

ALCM 

ft 


Length, confluence point to plane of 
skirt 

MI 

y 

CA PMAS 
GAMMA 

slugs 

rad 


Included mass 
Suspension line angle 

KX A1 

IXXA1 



Apparent mass tensor 

iyy ai 

IYYAI 



* Diagonal Elements 

IZZ A1 

IZZA1 


4 


IXX* 

IXX10 



Total parachute 

IYY* 

IYY10 



> Inertia Matrix 

IZZ* 

IZZ10 



Diagonal Elements 

Lp 

C M 

LCM 

ft 


Length, confluence point to plane of 
skirt 

c p 

LCP 

ft 


Length, plane of skirt to center of 
pressure 

LS 

LS 

ft 


Suspension line length 

LI 

LI 

ft 


LCM + LCP 

L2 

L2 

ft 


Riser length 

L3 

L3 

1 

ft 


SRB Center of Mass Location from 
nose 

L4 

L4 

ft 


SRB Center of pressure location 
from center of mass 

N 

M 

— 


Number of suspension lineb 

m C 

MC 

( 

slugs 


Canopy mass 

m L 

ML 

slugs 


Suspension lines mass 

m p 

MP 

slugs 


mj + m la 

m l 

Ml 

slugs 


rn + m T 
c L 

m la 

MIA 

slugs 


Apparent mass 

m 3 

M3 

slugs 


SRB mass 
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TABLE A 5 - LIST OF SYMBOLS FOR SUBROUTINE CHUTE 
(CONCLUDED) 


Quantity 

Mnemonic 

Units 

Description j 

0 

RHO 

slugs /ft 3 


Air density 

R o 

RO 

ft 


Skirt diameter 


SI 

ft 2 


Nominal area of parachute 

% 

S3 

«* 


Nominal area SRB 


XC 

slug ft 2 


^ Canopy moments of inertia 


ZC 

slug ft 2 

4 % 




XSL 

slug fr 


Suspension lines moments of inertia 


ZSL 

1 slug ft 2 

J 

) 


Coefficients for the tangent force are calculated as function of the angle of 
attack, a , using the polynomial form 

Cj » + C 2 a 2 + C3 a 3 + + C Q a B + C Q a 9 


Specifically for the parachute the normal force coefficient polynomial is of 
order three, the tangent force coefficient polynomial is of order five, and 
the moment coefficient polynomial is of order eight. 

The SRB normal force coefficient polynomial is of order eight, the tangent 
force coefficient polynomial is of order five, and the moment coefficient 
polynomial is of order nine. 


Angle of Attack -- The angle of attack is defined as the angle between the 
body axis of symmetry and the relative velocity vector. 



Sideslip Angle -- The side slip angle is defined for this problem to be the 
angle between the body fixed X axis and the projection of the relative velo- 
city vector on the body fixed X-Y plane. Thus, 
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0 ; 


Tan 


-1 



Va Xi 


Subroutine COEFTS is diagrammed in Figure A10 and listed in Figure All. 
Principal variables are listed in Table A6. 

Subroutine FORCES. Subroutine MOMENTS -- The subroutines FORCES 
and MOMENTS calculate the aerodynamic forces and total external (aero- 
dynamic and constraint) moments on the parachute and the payhxr . The 
dynamic pressure at the center of pressure of each body io calculated. 

Subroutine FORCES is diagrammed in Figure A12 and listed in Figure 
A 13, and its principal variables are listed in Table A7. 

Subroutine MOMENTS is diagrammed in Figure A 14 and listed in Figure 
A15, and its principal variables listed in Table A8. 

Subroutine DIRCOS. Subroutine DBDT -- Subroutines DIRCOS and Sub- 
routine DBDT calculate ancl manipulate the matrices of direction cosines 
describing the orientations of the reference frame, parachute, riser, and 
payload with respect to the earth. DIRCOS calculates the immediate direc- 
tion cosines matrices as functions of the Euler angles at each integration 
step. 

For resolution of the riser force (the constraint force) into the parachute 
and payload reference frames direction*, d i ction cosines matr s are 
formed describing the orientations of - fixed axis system . ith 

respect to the parachute and the payload oudy a..:ed axes systems. 

Subroutine DIRCOS is diagrammed in Figure A16 and listed in Figure A 17, 
and its principal variables defined in Table A9. 

Subroutines DBDT is d*agrammed in Figure Aid and listed in Figure A 19, 
and its principal variables defined in Table A 10. 

Subroutine PRECOR -- Subroutine PRECOR integrates the equations of 
motion using a Runge Kutta initialization and a predictor-corrector integra- 
tion algorithm (Ref. 12). 

The Runge Kutta method establishes values for the state vector at time 
zero and at time equal to one integration step size. Using these *wo initial 
points the state vector is updated in the predictor mode (mode = 5) and 
time is increased one integration step size. The corrector mode (mode - 
6) refines the prediction made when mode = 5. Completion of the correc- 
tions returns control to the main program for calculation of everything 
associated with the newly calculated state vector. 
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Figure A10. Subroutine COEFTS Flow Diagram 
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SUBROUTINE COFFTS 4RCI 

C SUBROUTINE COFFTS CALCULATES THF AERODYNAMIC COEFFICIENTS OF TMF PARA- C 

c chute and thc srb as functions of alti iude*velocity* and anglf of attack c 

COMMON/ AAA/ ACT 4 9 > • AfN 4 9 1 * ACM! 9 1 *RCT C 9 1 »RCN I 9 1 • RC* ( 9 ) 

COMMON/AAN/ Y433) 

COMMON/AAD/Rf 9»3»9>»ft546>»T<3»6> 
r COMMON/AAF/CNX »CTI »CN9*CT9 tALPHAl •ALPHA9*RFTA1 •RFTA3»GAMMA 

COMMON /AAFF/ CM 1 *CM3 *ALCM 

C0MM0N/AAG/L2#L20OT.L2OD0T#LCN t LCNDOT#LCNDnT 
COMMON/ A AH/C 1 #C3*F2*L3*RAD#L1 *LA*CF l*CF3*Sl*K‘* 

COMMON /AAP/YW I NO* VW I ND # VGUST • W I Gil 
REAL tl #L3*L4 

REAL L2*L2OOT#L200OT#LCM#LCM0OT»LCM0DT 
C VELOCITIES SQUARED OF THF CCS OF »OOIFS 1 AND 9 
C INERTIAL VELOCITIES SOMARFD OF THC CM»S OF BODIES 1 AN* 3 
ci * Ym*Yn»*y(?)#yj?)+Yfi)*Yn) 

C3 * Yf 10l*vn0)4.Yltl |#YM 1 14.YU2) *YC1?| 

C VELOCITIES SQUARED AT THF CP»S OF BODIES 1 AND 3 WRT T HF AIR *ASS 

CF1 * 4Y4 l >-VI5>*ILl-LCMl-WIGU*BIl*l#in**?<MY4?>*Y|4l*4Ll-LCM> 

1 ~WIGU*B4 1 •?*! ) I **?+! Yf 3)— WlGU*P4 1 1 3 »1 ) ) **2 

CF3 * CYf 10l+VtlA)*LA-NlG4^Pf^,l*ln**2+|YCn)-Yf 13)*L4-u*IGu*P< 

1 3*?.m**?+«YU?)-WIGU*o(3*3,in**? 

C NORMAL COMPONFNTS OF VELOCITY AT THE CP OF POO IFF 1 AND 3 
C WITH RFSPFCT TO TMF AIR MASS 

CNORF1 * $ORT|ARSlCFl-IY43»-WlGU*R| 1*3*1))**?)) 

CNORF* * SQRTCABS«CF3-IYCl7)-WIGU*BC‘*,3*lll#*?n 
C ANOLFS OF ATTACK 

ALPHA 1 * ATANfCNORr 1/|y(3)-WIGU*F{l*3 # n + l.-F-l<oj 

IF 4 V 4 1 )—Y 4 5)* ( Ll—LCM )-wlGU*R4 ItUlUGT. o.'M ALPHAI * — alpha i 

ALPHA 3 * AT AN 4CNORF3/ 4 Y4 12 )~W!GU4D4 3*3*1 1 + 1 *OF“l A I ) 

IF 4Y|10)+y4 l4)*L4-WlGU*«43*ltl)»GT*0*0) ALPHA 3 * -ALPHA3 
C ANGLE BETWEEN THE X AXIS AND THE PROJECTION OF C ON THF X-Y PLANF* 

C BODIES 1 AND 3 

BETA1 * ATAN4 4YC?)+Y4A)*4Ll-LCM)“WlGU*B41*?*t M/4Y4 U-Y 4 5 I* 4L 1-LC 
1 M 1 *W T OU*° 4 1 • 1 • 1 ’^nOF-TA) ) 

BFT A3 « ATANf 4Y41 1 )-Y 4 1 3 >*L4-WlGU*Bf 3 * 2 * 1 1 I / 4 V 4 10)+Y( 14 )*LA~WIGU* 
l«43 t l#n^l*OF*lAM 
C PARACHUTF NORMAL FOPCF COFFFICIENT 

'•fcil = AC M ( 1 )*ALPHA1 ♦Arw{?)#ALPHAl*«>+ACM(3)*ALPHAl*#^ 

^ fpn *»OD»Mt FOPCF rOFFFfr?FNT 

CN3 * BCNCl )*AlPHA34-PCN(?)*ALPHA9*#?^nc^f^I *ALPHA3**^ 

1 +BCN 4 4 ) *AlPHA 3*#4+RCN C f ) *ALPH A3**5+BCN C 6 ) *ALPHA 3**6 

2 ♦BCN 4 7 ) 9ALPHA 3»*7+Bt N f 8 ) *ALPHA****8 
r PARACHUTF TANGFNT FORCF COEFFICIENT 

CTI * ACT 4 1 I ♦ACT 4 ? ) *ALPhA 1+ACT 43) * ALPHA t #*?♦ ACT C 4 } *ALPHA 1»»'» 

1 +ACT(5 )*ALPHA1##4+ACT4 6)*ALPHA1#*«; 

r ff QpCP COFFFTrfFMT 

CT3 * BCT( 1 ) ♦BCT 4 ? ) *ALPHA3*RCT ( 3 ) *AL PH A i •♦/♦BCT 44) ♦ALPNA3**'* 

1 4BrT(S)*ALPHA3* # 44flCT4A)*ALPMA3**S 

? +RCTC7)*ALPHA3**6«‘RCT(8>*ALPHA3**7+RCT49)*ALPHA3**8 

= AC^nJ^AtPHAl 
rM-a * °CM( 1)*ALPHA3 

00 to 

C PARACHUTF MOMENT COFFFICIFNT 

CM1 » CM14ACM4I JMLPHAl**! 

C PAYLOAD MOMENT COFFFICIFNT 

]0 CM3 « )*ALPHA9**I 

0ff?(f9A? 

r mo 


Figure All. Subroutine COEFTS Source Listing 
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TABLE A6 - LIST OF SYMBOLS FOR SUBROUTINE COEFTS 


Quantity 

Mnemonic 

Units 

Description 


ACM 

--- 

Constants in polynomials 



ACN 

— 

for parachute aerodynamic 



ACT 

— 

coefficient 


a l 

ALPHA 1 

rad 

Parachute angle of attack 


*3 

ALPHA 3 

rad 

SRB angle of attack 



BCM 

— 

| Constants in polynomials 



BCN 

— 

I for SRB aerodynamic 



BCT 

— 

coefficients 


h 

BETA 1 

rad 

Parachute sideslip angle 


% 

BETA 3 

rad 

SRB sideslip angle 


Va l 2 

CF1 

(ft/sec 2 ) 

Velocities squared of the 

parachute 




CP's WRT 


Va 3 2 

CF3 

(ft/sec 2 ) 

the moving air mass 

SRB 

C M t 

CM1 

--- 

. 

parachute 

1 



Moment coefficients | 


° M 3 

CM3 



SRB 

C N- 

CN1 

--- 


parachute 

1 



Normal force coefficients 


C N 

N 3 

CN3 



SRB 

C T« 

CT1 

— 


parachute 

*1 



Tangent force coefficients 


Ct 

*3 

CT3 



SRB 

c i 

Cl 

(ft/sec) 2 

Inertial velocities 

1 parachute 

C 3 

C3 

<ft/sec) 2 

Squared 

| SRB 
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Figure A12. Subroutine FORCES Flow Diagram 
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subroutine forces 

c subroutine forces calculates the aerodynamic forces on the parachute and srr c 

C AS FUNCTIONS OF THE VFLOCI.Y OF THE CP RFLATIVF TO THE AIR AND THE ANGLF OF C 
C ATT AC* C 

COMNON/AAF/CNl »CTX • CN3 »CT3» ALPHA 1 »ALPHA3,BFTA1 • RET A3 t GAMMA 
COMMON/ A AH/C l *C3»F2»L3»RA0»L1 »LA*CF 1 »CF3» Si *S3 
COMMON / A AL/F 1X*F1Y»F1Z»F3X»E3Y*F3Z 

C0MM0N/AAR/KLS»KR»MC*ML»MI»MlA»MPtD0»R0»RH0»L?0*LCM0*LS0»M*M3 
RFAL LI *L3*L4 

REAL KLS*KR*MC*ML.'M .Ml A. mP»L?0»LCM0.LS0.M,M3 
C DYNAMIC PRESSURE, BODIES 1 AND 3 
OSI • O.S*CF1*SI»RHO 

OS3 » O.S»CF3*S3*RHO 

C AEROOTNAMIC FORCES IN X* Y. ANO Z BODY EtxEO AXIS DIRECTIONS* BODIES I AND 3 

MX * ♦CNl*OS]*rOSt«FTAl ) 

FIY « ♦CN1»OS1#SIm(«»FT* 1 1 

FIT « -CT1*0S1 

F3X « ♦CN3*OS3*COSI»FTA31 

F3Y • ♦CN3«OS3*SlNI°FTA3l 

F3? • -CT3*OS3 

RETURN 

FND 


Figure A13. Subroutine FORCES Source Listing 
TABLE A 7 - LIST OF SYMBOLS FOR SUBROUTINE FORCES 


Quantity 

Mnemonic 

Units 

Description 

F 1X 

FIX 

lb 

\ parachute aerodynamic 

F 1Y 

FIY 

lb 

\ forces in XYZ 

F IZ 

FIZ 

lb 

j body fixed axes directions 

F 3x 

F3X 

lb 

^ SRB aerodynamic forces 

F 3Y 

F3Y 

lb 

> in XYZ body fixed 

F 3Z 

F3Z 

lb 

) directions 

q l St >! 

QS1 

lb /ft 2 

1 v 2 c 
2 pV al S ° 1 

q3 s o 3 

QS3 

lb/ft 2 

1 v 2 c 

2 p V a3 S o 3 
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Figure A 14. Subroutine Moments Flow Diagram 
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Figure A15. Subroutine Moments Source Listing 
TABLE A8 - LIST OF SYMBOLS FOR SUBROUTINE MOMENTS 


Quantity 

Mnemonic 

Units 

Description 

M 1X 

MIX 

ft-lb 


total external moments 

“lY 

M1Y 

ft-lb 


, about XYZ parachute 

M 1Z 

M1Z 

ft-lb 


body fixed axes 

M 3X 

M3X 

ft-lb 

» 

total external moments 

M 3Y 

M3Y 

ft-lb 


> about XYZ SRB 

M 3Z 

M3Z 

ft-lb 

4 

body fixed axes 
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Figure A16. Subroutine DIRCOS Flow Diagram 
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SUBROUTINE DlltCOS 

C<HHMHHHHWmhhhmmmm. M «n»AN W— ' a m aara a aa r — AAAAAAAAA f AAAi W MA M AA AA AAAl 

c SUBROUTINE OIRCOS (DIRECTION COSINES) FORMS_THE JNMEPlATEOIRECTipN, 

C COSINES MATRICES FOR THE PARACHUTE. fHE RI SER S, AND THE PAYLOAD. AMO 
C FORMS THE RELATIONS NEEDED TO R£ SOLVE tHf FORCE IN THE RISER INTO THE BODY 
C FIXED COORDINATE SYSTEMS OF THE PARACHUTE AND THE PAYLOAD _ _ . .. 

C THE DIRECTION COSINE MATICES ARC NOTED AS FOLLOMS BU.J.KI MMERt I IS 
C THE PARTICULAR REFERENCE FRAME. J IS THE RON NUMBER. AND R IS THE 
C COLUMN NUMBER FOR THE ELEMENTS IN THE 3X3 MATRIX. __J.. 


COMMON/AAB/ Y(33) 

COMMON/AAO/Bt 3.3.3). BS(S).TO.S) 
X ■ 1 

DO 20 IJ * 1.3.2 
I • IJ 

IF (IJ «EO*S) I > 2 
X - Xa6 

DO 10 J » 2.6.2 
TII.J-1) • SINCVtX) ) 

TII.J) » COSIYIX) I 

10 X » XAl 

20 CONTINUE 

00 30 U » 1.3.2 

1 • IJ 

IF (IJ .EO.S) I - 2 


Bit. 1.1) 
Bit. 1.2) 


T ( I »6)*T ( I »A) 
TII»5)*TI1»A) 


B( 1.1.3) — TII.3) 


B( I .2.1 ) 
BII.2.2) 
BII.2.3) 
B( I .3.1 ) 
Bll.3.2) 
30 Bl I .3.3) 
X 

DO AO I 
00 AO J 
BS(K) 

AO X 

RETURN 

ENO 


TII.l)*TII.3)*TtI.6)-T(I.2)*T(I»5) 

T ( I .5) *TI I. 3) *T( I *1)ATII .6|*T( 1.2) 

Td.A)AT(I.l) 

T(I. 6) *T II. 3 )*T (1.2) AT (I* 5) *T (1.1) 

T(I.2)AT( I.3)*T( I*S)-Tlt.l)*t(Ii6) 

Td.A)AT(t.F) 

1 

1.3.2 

B ( I .J.l ) *8(2.3. 1 )aB( I.J.2)*B(2»3.2)AB(I . J.3 )*B( 2.3.3) 
XAl 


Figure A 17 . Subroutine DIRCOS Source Listing 
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TABLE A9 - LIST OF SYMBOLS FOR SUBROUTINE DIRCOS 


Quantity 

Mnemonic 

Units 

Description 





Direction cosine matrix 
elements i, k = 1, 2, 3 

4 

B(J, I, K) 



j = 1 parachute 
j = 2 riser 
j = 3 SRB 

B S1 

BS(1) 


/ 

used for rotating a vector in 
Earth coordinates to one in 
j coordinate system 

proportion of Fg projected 

B S2 

BS(2) 

— 


on X, Y, Z parachute 

B S3 

BS(3) 

— 


1 body-fixed axes 

B S4 

BS(4) 

— - 


proportion of F 2 projected 

B S5 

BS(5) 

— 


on X, Y, Z SRB body 

B S6 

BS(6) 

— 

l 

fixed axes 
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Figure A18. Subroutine DBDT Flow Diagram 


141 
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♦ ? 


AC ( 

T+f )* 

Y ( *> ) 

#n 

( 1 ^ 

»n 

• V 

4) *0(1 

*? 



DPT! »0*‘ 
rwr> 


Figure A19, Subroutine DBDT Source Listing 


TABLE A10 - LIST OF SYMBOLS FOR SUBROUTINE DBDT 


Quantity 

Mnemonic 

Units 

Description 

*■ * 

AA (array) 

• * ■“ 

Array containing elements 
of the first columns of the 
time derivatives matrices 
of bodies 1, 2, and 3 

— 

AB (array) 

— 

Second column elements 

— 

AC (array) 

— 

Third column elements 
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Subroutine PRECOR is diagrammed in Figure A20 and listed in Figure 
A21, and its principal variables are listed in Table All. 

Subroutine WIND -- Subroutine WIND calculates at each integration step 
the value of the 5% risk wind speed profile as a function of the altitude. 

The wind velocity vector is assumed to be aligned with the earth-fixed 
reference frame X axis. 

Subroutine WIND is diagrammed in Figure A22 and listed in Figure A23, 
and its principal variables are listed in Table A12. 

Subroutine GUST -- Subroutine GUST computes a step change in the air 
mass velocity vector according to a 5% risk gust envelope related to the 
5% risk wind profile. The step changes are calculated at a frequency of 
four per minute of simulation time and are both sign and magnitude modified 
by a random function. 

Subroutine GUST is dia£ ammed in Figure A24 and listed in Figure A25, and 
its principal variables are listed in Table A 13. 

Subroutine ELASTIC - - When the elasticity option is employed, subroutine 
ELASTIC is called at two-second intervals to determine the first and 
second time derivatives of the lengths of the elastic elements, the riser, 
and the suspension lines. The method employs a central difference method 
on an averaged length. 

Subroutine ELASTIC is diagrammed in Figure A26 and listed in Figure 
A27, and its principal variables are listed in Table A16. 

Subroutine PRINT -- Subroutine PRINT controls the line printer operation 
and loads plotting storage arrays. Ten groups of data are printed on each 
page. This is adjusted by changing the line output counter (LOC). When the 
number of groups printed equals LOC, a heading is printed at the top of the 
next page and the LOC is set to zero. 

Corresponding to each output group, the values for altitude, range, angles 
of attack, pitch angles (0i), riser force, riser length, center of parachute 
mass, and the air mass velocity are loaded into arrays for use in plotting. 

Subroutine PRINT is diagrammed in Figure A28 and listed in Figure A29. 

Its principal variables are listed in Table A 15. 

Subroutine CONST --S ubroutine CONST calculates a group of variable com- 
binations used in the differential equations subroutine DIFEQN that result 
from the method of coupling of the parachute and payload. Generally, these 
are the accelerations of the confluence point and attach point in components 
parallel to the earth fixed axis system. 

Subroutine CONST is diagrammed in Figure A30 and listed in Figure A31. 
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aim^romtirf prfcor <*,f»m> 


<i»«* 0»IT!MF PR p rOR CPRrOTrToR - CORRFF TOR) IRTFORATF* TMF SYSTFH OF 
niPFFRFRTl^L FCUATICN^* IT* FFATURF4 TNCLUOF A Ri#IC%A«MCUTTA tN*T I ALI2ATT0R* 


C»**w 


1 


* 


rORROR/AAR/ YI33I 

CORROR/AAJ/’-OOF 

rORW>H/AAT/TT M P 

0TRFR5I0N nl33>*AOn3>*AFt33|*FI33|*PVC33>*CYl33)*0C33l*SC33l 
ROOF * ROOF ♦ 1 
TF (ROOF *OT* 61 ROOF » 6 
rORTfRUF 


M * .01*62* 

tP (TfRF*LT*T**l h * 0*007*1?* 


?r |TfRF # r.r*?*.Pi m * o*oai?* 
OOTO n*4»«t 0*14*19) » ROOF 
TfRF* * TfRF 

"O A T * It* 

*U> * Ylt) 

Oil) * 0*0 

""m * vr r > 

ocn * Frn*H 

Y( n * Ml) ♦ o.A*h*FU) 

TIMP « Tt^F* +0*3»H 


tl 

a *n f *l*R 

Yf f > * t(»U^,S#H#F(TI 

ocn * ocn+7.o*H*Fm 

AO 0(|| * P(f| 

OOTO 13 

• A *0 4r» t*f*« 

vfM « A(tl+K«F(li 

40 om * 0!Ti4?*o#H»rn) 

TfRF « TfRF* ♦ H 

wvron 

9 oo 70 1*1 *R 

yo ycti ■ *cn*coc n+Fm*w>/6.o 

n r * 1 

Oct'iq*' 


1A *0 IT T* 1*N 

r PRFnfrroo rooATTO* 

PYf n • AOfl )*?.0*N*Ff tt 

*m * fui 
AF(n » ycij 
17 V(f ) * PYCtl 

TfRF « TJRF^H 
J * n 

to r • i 
RFTtiR” 

1 P J * J4l 

OO TO f« 1 »W 
C rORRFCTOR FOUATION 


CY( f I • AF< n*0*9*H*fO(! )♦?(!») 

20 Ytn • cyu) 

TF U*FO*?| r>OTO 3? 

RFTtlRR 

32 DO 100 I *1 #N 
AOffl • AFC 1 1 
100 Ml) » YIU 
T1RF5 ■ TTRF 
MOOF • A 
RFTfJRN 
FRO 


Figure A21 


Subroutine PRECOR Source Listing 
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TABLE All - LIST OF SYMBOLS FOR SUBROUTINE PRECOR 


Quantity 

Mnemonic 

Units 

Description 

— 

ATHl) 

— 

Y(J)| t-H 

— 

AE(I) 

— 

Y(I)|t 



CY(I) 

— 

corrected value Y(I) j t+H 

— 

DU) 

— 


— 

F(I) 

— 

s * H V h 

At 

H 

sec 

stepsize 

— 

MODE 

- — * ■ 

MODE = 4, Runge Kutta initialization 
MODE = 5, Predict MODE = 6 
correct 

— 

N 

— 

number of equations 

— 

PY(J) 

— 

predicted value Y(I)| ^ 

t 

TIME 

sec 

time 

— 

Y(I) 

— 

state vector 
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Figure A22. Subroutine WIND Flow Diagram 
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SUPROUTINF WINO 

C SUBROUTINE WIND CALCULATES A WINO INPUT TO THE VELOCITY OP THE AIR HASS C 

COHNON/AAB/ Y(3TI 
COMMON/AAP/YWIND*VWlNDtV6UST.WI6U 
IF (YI30J.LT.-49S.0) 6OT0 10 
VWINO ■ 69«0*C— YI301/A95»01**0#21 
GOTO 20 

10 VWTNO - 69.0 

?0 CONTINUE 
BFTURN 
END 


Figure A23. Subroutine WIND Source Listing 


TABLE A 12 - LIST OF SYMBOLS FOR SUBROUTINE WIND 



Mnemonic 

Units 

Description 

v 

wind 

VWIND 

ft/sec 

mean wind speed 
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Figure A24. Subroutine GUST Flow Diagram 
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SUBROUTINE GUST ' 

C SUBROUTINE GUST CALCULATES A GUST INRUT TO THE VELOCITY OR THE OMYONOT OM» 

C MASS AS A FUNCTION OF ALTITUDE A NO A RANDOM MAGNITUDE MODIFIER 

COMMON/AAR/ YC33J ' 

COMMON/ A AP/YW I NO* VW I NO • VGUST • W I GU 
20 IFIYC30I.LT. -990.01 GOTO 30 

VGUST • 19.8 v. v 

GOTO 50 

30 IF ( YC 30 I *LT.-3280.0)GOTO *0 

VGUST * 10. 8 /2290. 0# C -Y ( 301 -990.CT) +H.9 

GOTO 50 

AO VGUST • 29.7 
30 IEIKr.EO.Ot SEED • 7.0 

KK -1 *\< 

CALL RANDU (YFL.SFEol 

C YFL IS IN THF RANGE -1.0 to 4-1*0 

VGUST ■ YFL*VGUST 

RETURN 

END 


Figure A25. Subroutine GUST Source Listing 


TABLE A13 - LIST OF SYMBOLS FOR SUBROUTINE GUST 


Quantity 

Mnemonic 

Units 

Description 

V GUST 

VGUST 

YFL 

ft/sec 

gust velocity 

random modifier 
in range -1 s YFL * 1 
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SUBROUTINE ELASTIC 

C SUBROUTINE ELASTIC CALCULATES THE TINE RATE Of CHANGE OF THE LENGTH OF TH C 

C ELASTIC RISERS AND SUSPENSION LINES AS WELL AS THE TIME RATE OF CHANGE OF C 

C THE RELATIVE VELOCITIES OF TSCH END OF THE ELASTIC ELEMENTS C 

COMMON/AAG/L2 .L2DOT.L2DDOT .LCM. LCMDOT .LCMDOT 
COMMON/AAK/AL1.AL2.AO.AL4 
COMMON/AAO/OLDTINE»TELAST*ET!ME 
C OMM O N /AAT/Tt HE 

COMMON /XOROS /SUHMA 1 »SUMMA2 . TOTAL .AVER A 1 * AVER A2 .OVER A 1 . OVER A2 » OT 
REAL L2.L2DOT.L2DOOT.LCM.LCMDOT .LCMDOT 
C ELASTICITY CALCULATIONS 

L20OT. (AVERAI-OVERAlt/DT 

LCMOOT* IAVERA2-OVERA2I/DT 

L2000T- (AL2-L2DOT)/OT 

LCMDOOT» (ALA-LCMOOTI/OT 

AL1 • L2 

AL2 • L2DOT 

ALT ■ LCM 

ALA • LCMDOT 

OVERA1* AVERA1 

OVERA2* AVERA? 

RETURN 

END 


Figure A27. Subroutine ELASTIC Source Listing 
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TABLE A 14 - LIST OF SYMBOLS FOR SUBROUTINE ELASTIC 



Mnemonic 

Units 

Description 

— 

AL1 

ft 

Last calculated value Lg 

— 

AL2 

ft /sec 

• 

Last calculated value Lg 

— 

AL3 

ft 

Last calculated value L n 

C M 

— 

AL4 

ft/sec 

Last calculated value L„ 

^M 

— 

A VERA 1 

ft 

Average value of Lo during interval 
from t - DT to t - DT/2 

* * ** 

AVERA2 

ft 

Average value of h r during interval 

M 

from t - DT/2 to t 

2 At 

DT 

sec 

Averaging interval 


LCM 

ft 

Length from confluence point to para- 
chute center of mass 

l c 

LCMDOT 

ft/sec 

■4 L_ at » - DT/2 

® C M 

L 2 

L2 

ft 

Length of riser 

l 2 

L2DOT 

ft/sec 

d L 

dt l C m at t - DT/2 

— 

OVERA1 

ft 

Average value of L 2 during the inter- 
val from t - DT to DT/2 

~ ™ * 

O VERA 2 

ft 

Average value of h r during the 

C M 

interval from t - DT to t - DT/2 

t 

TIME 

sec 

time 
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Figure A28. Subroutine PRINT Flow Diagram 













subroutine print 

C SUBROUTINE PRINT HANDLES ALL THE PRINT OUTPUT FUNCTINS OP THE PROGRAM 
COMMON/AAP/ VI 33) 

COMMON/AAC/D! 30) 

COMMON/AAD/B * 3*3*3) *AS(6t.T(3*6l 

C0MM0N/AAF/CN1»CT1#CN3*CT3»ALPHA1 *ALPHA3*BETA1. BETAS .GAMMA 

C0MM0N/AA6/L 2 .L2D0T »L2DD0T .LCM.LCMDOT *LCNODT 

COMMON/AAH/C 1 »C3»F2*L3»RAD»L1*L4»CF1»CF3* 51*53 

COMMON/AAP/VWINO.VWIND.VGUST,KIGU 

COMMON/AAT/TIME 

COMMON/AAW/NORS 

COMMON/PL TR/XXI 402 ) »TI(E1 (402) *THE3( 402) .API (402) .AP3I402) .ALT (402) 
l.RN6(402).FOR(402).RL(402).CL(402).W6(402) 

REAL L2.LCM.L3.L1.L4 

REAL L2DOT.L2nl>OT.LrMOOT.LCMODT 

IF (TIMF.FQ.O.O) LOr * 10 

VI 31 ) « -Ml.3.1)*LCM-R(7.3.1)*L2-BI3.3.1)*L34Y(28) 

VI 33 ) « -Btl,3.3)4LCM-R(2.3.3)*L2-BI3.3.3)*LS4Vt30) 

VI 32) * -B(1.3.2)*LCN-e(2.3.2)*L2-BI3.S,2>*L3 +YI29) 

CALL T00E6 
FI « S0RTIC1) 

F3 » SORT (C3) 

ALP1 ■ ALPHA 1*RAD 

ALPHAS » ALPHA3FRA0 

IF (LOC .FO. 10) SOTO 300 
GOTO 330 
300 IOC « 0 

IF (TTMF.FO.0,0) LOr * 10 
WRITE (6.530) 

330 FORMAT (1H1) 

320 WRITF (6*300) 


XX (NOBS) * 

T IMF 

TWF1 (NOBS) 

* 

V ( B ) 

THF3(N0RS) 

a 

V ( 1 7 ) 

API 

(NOBS) 

a 

ALP1 

AP3 

(NOBS) 

a 

ALPHA3 

ALT 

(NOBS) 

a 

•Yf 30) 

RNG 

(NOBS) 

a 

Y(28) 

FOR 

(NOBS) 

a 

F? 

RL 

(NOBS) 

a 

1? 

r t 

(NOBS) 

a 

LCM 

WG 

(NOBS) 

a 

WI GU 


133 CONTI NUF 
CALL TORAO 
ALP1 « ALP1/RA0 
ALPHA3 ■ ALPHA3/RA0 

400 FORMAT! 10X.SHVWIND.F8.3. 5X.SHVGUST.PB. 3./) 

300 FORMAT (23X*9HPARACHUTE»43X*3HSRB/2X»4HT!Me*4X*4HPStl* 

16X*2HU1*7X»2HP1*TX*3HXE1 *6X*2HC1*TX *4HPSI3»3X»2HU3»TX*2HP3 
2.7X.3HXE3.5X *6HXE3D0T.4X»2HC3.7X.3HLCM.6X»2HP2/9X»6HTHETA1 
3*3X»?HV1 *7X* 2H01 »7X»3HVE1 *14X*6HTHETA3»4X*7HV3*7X*2H03*7X* 
43HYF7.3X,6HYF3DOT.13X.7ML2/10X.4HPHI1.6X.2HWI.TX.2HR1*TX»3HZE1. 
34X.6HAL»HA1.SX.4HPHI3.SX.2HW3.7X.2HR3.7X.3MZE3.SX.6HZE3D0T. 
63X.6HALPHA3. / ) , 

310 F0RMATIP7.2. 12F9»3»F12.S/TX*4P9.3*9X»SP9»3*9X»P9*3/7X»11F9»3 ) 
RETURN 
FNO 

330 WRITE 16*310) TIME.V19) .VI 1) *V(4) .VI31 ) * El .VI 18 ) *V( 10) *V( 13 ) . 

1VI28).0C78)*E3.LCM*F?.VI8)*VI2).VI3)»VI32), VI 17) *VI 11 ) *VI 14 ) . 

2VI29) *0129) *L2* V(7) »V(3) »V(6) .VI33) .ALP1 .VI 16 ) »VI 12) .VI 19) . 

3V( 30) *0(30) .ALPHAS 
WRITE (6*400) VWIND.VGUST 
LOC • L0C-.1 
C LOAO PLOT ARRAYS 
C UP TO 400 POINTS PPR CtlRvF 

IFtNOBS.0E.400) QOTO 139 
NOBS - NOBS+l 


Figure A29. Subroutine PRINT Source Listing 
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TABLE A15 - LIST OF SYMBOLS FOR SUBROUTINE PRINT 


Quantity 

Mnemonic 

Units 

Description 

*1 

ALP1 

deg 

parachute angle of attack 

h 

ALT 

ft 

plotting storage array, altitude 


API 

deg 

plotting storage array parachute 
angle of attack 

a 3 

AP3 

deg 

plotting storage array SRB angle of 
attack 


CL 

ft 

plotting storr . ;e array, Lp 

M 

VcT 

El 

ft/sec 

speed, parachute center of mass 

Vs 

E3 

ft/ sec 

speed, SRB center of mass 

F 2 

FOR 

lb 

plotting storage array, riser force 

— 

LOC 

— 

line output count 

— 

NOBS 

— 

number of points in each curve 

4 

RL 

ft 

plotting storage array, riser length 

XE 3 

RNG 

ft 

plotting storage array, range 

*1 

THE1 

deg 

plotting storage array, 0^ 

*3 

THE3 

deg 

plotting storage array, 0g 

— 

WG = V/IGU 

ft/ sec 

air mass velocity vector 

t 

XX 

sec 

plotting storage array, time 


158 











c 



Figui e A30. Subroutine CONST ^low Diagram 
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SUP ROUT INF CONST 

C SUBROUT INF CONST CALCUL ATFS CONSTANTS USFH IN TMF 0 IFFFRFNTI AL fCUaTIONS 
COMMON/AAr,/U?*L?OOT»L?OOOT»LC^*LCVDOT »LCMD n T 
COMMON/AAH/C1 ♦C3»F?*l3«RAr>*Ll ♦L4*CFl«CFT*Sl 
COMMON/AAP/YI S3) 

COMMON /AAC/O ( SO ) 

COMMON/ A AD/B (3»3*3)*PS(6)*T<3»6) 

COMMON/AAN/A6»A7*A8»A9»A10*AU«A1?»A1T.A1A,A1 A,A1S* A17,a: fl»AI0,A?n 

COMMON/ AANN/A?1 

RFAL L2«L2r>0T,L20D0T.Ll .L4 


RFAL L3»LCMtLCMr>0T*LCMf>r>T 

A6 


(D( 10l-D(l4)»H)*(M(3*l»3> > 

A7 

x 

<D< 11)*0(1‘»>*L3)*“(3»?»3> 

AB 

* 

(Oil OtB1*LCM+YIS)*LCMnOTJ»B( 1»’ »? ) 

A9 

X 

C 0 ( 2 1 -D ( A ) *LCM-Y ( A ) ♦LCMHOT >*"(].:>,•») 

A10 

s 

(0( l?)*F>n*3*7l ) 

All 

X 

(D( 10l-ru H)*L3)*»(3»1*?> 

A12 

« 

(D( 11»+D(13)*L3)*«<3«'’ 21 

A 13 

X 

D(12)*B(3»3.?i 

A1A 

X 

(Oil )-*>0(5)»LCM+Y(51«LC.(OOT>»R( 1 • 1 • ? 1 

Al* 

X 

( D( 2 ) “Ot 4 ) *LCM—Y ( A ) ALCMDOT )*B(1#2»2I 

A16 

X 

(Dnoi-f>CA>*L7)*"(?tl »1 ) 

A17 

e 

(f>( U>+0( l’I*L3)*“n»2*l) 

Al 0 

X 

D(12)*H!3»S.1> 

A*9 

X 

( D< 1 »♦!>( 9 )*».CM+Y( 5 l»LCMOOT I »P ( 1 . 1 , 1 1 

A20 

X 

(0(2 )“D(4 ) *LC M ~Y ( A )*LCMDOT ) *R ( 1 » 2 * 1 ) 

A?1 

RFTURN 

FND 

X 

DCM+LCMOm 


Figure A31. Subroutine CONST Source Listing 
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Subroutine INVELO -- Subroutine INVELO initializes the inertial components 
of velocity in the body fixed axis systems at time zero for the initial orien- 
tations and vertical rate of descent as read in the input data deck. 

Subroutine INVELO is diagrammed in Figure A32 and listed in Figure A33. 


Auxiliary Subroutines 

Auxiliary trigonometric functions SEC provided. SEC is listed as Figure 
A34. 

Subroutine SRBIN calculates the SRB inertial differences as used in Equa- 
tion (8A). Subroutine SRBIN is listed as Figure A35. 

Subroutine TORAD converts angles and angular velocities to radians and 
radians per second. Subroutine TORAD is listed as Figure A36. 

Subroutine TODEG converts angles and angular velocities to degrees and 
degrees per second. Subroutine TODEG is listed as Figure A37. 

Subroutine RANDU calculates a random number in the range -1 to +1. Sub- 
routine RANDU is listed as Figure A38. 


Linearization Subroutines 

F /e subroutines make up the package to linearize and find the eigenvalues 
for a set of nonlinear differential equations. 

Subroutine DERIVE calculates the first partial derivative matrix. Sub- 
routine DERIVE is listed as Figure A39. 

Subroutine EIGEN is called from subroutine DERIVE and performs the 
control, storage, and output functions for the eigenvalue calculation pro- 
cess. Subroutine EIGEN is listed in Figure A40. 

Subroutine HESSEN is called from subroutine EIGEN and manipulates the 
matrix of first partial derivatives into the upper Hessenberg form. Sub- 
routine HESSEN is listed in Figure A41. 

Subroutine QRCALL is called from subroutine EIGEN and hence calls sub- 
routine QR. QRCALL is a double iterative eigenvalue approximation 
method using a quotient reduction sche 'e provided by QR. Subroutine 
QRCALL is listed ' Figure A42. and subroutine QR is listed as Figure 
A43. 
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SUBROUTINE IN VELO (HOOT I 

C SUBROUTINE INVELO CALCULATES THE INITIAL VELOCITIES OF THE CHUTE AND SR* 
COMMON/AAR/ Y(33) 

CONNON/AAD/A (3*3*3) »°S( 6 1 »T( A*6 ) 

call oircos 

C INITIAL INERTIAL VELOCITIES 
C Ul « Y(l» VI * Y(?) Wl « Y(3l 

C U3 * YllO) V3 * Y( 1 1 | W3 * Y(12J 

DO 30 1*1*3 

Y(II * HDOT*B( 1*1*3) 

Yf 1*91 « HOOT*B(3*I*3» 

30 CONTINUE 
RETURN 
END 


Figure A33. Subroutine INVELO Source Listing 


REAL FUNCTION SEC (X) 

SEC » l*0/(COSfX)+l«QE-lM 

RETURN 

END 


Figure A34. Function SEC Source Listing 
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Cl loom IT t fiF CRnJM 

C II i« ROllT IMt SRRIN CALCULATES THF INFRTIAL CHARACTFRI STICS OP THP SR® 
PDMMDm/A AP/T XXI * T YY 1 ♦ 1 7 7 1 »IXX1»JYY^»J77^ 
rnMMO»i/AAOO/ IXZ1 »I YX1 ,T?Y1 ,1X7* ♦! YX3 *IZY3 
RPAL tXXI »!YY1» I 771 • I XXI,! YY3. 1 Z7* 

RPAL l XZl » l YX1 »I 2Y1 # IXZ3 »I YX3 ♦ IZY3 

r INFRTIAL CHARACTERISTICS, SRP 
IZY3 « IZZ3-IYY3 

TYX3 « IYY1-IXX3 

|T73 « IXX3-I77* 

RPTMRN 

FND 


Figure A3 5. Subroutine SRBIN Source Listing 


SU«*RO»»TINF TORAD 

C SUBROUTINE TORAO CONVERTS ANGLES AND ANGULAR VFLOCITIES TO RADIANS 
CONMON/AAP/Y(33) 

COMNON/AAH/C1 »C3»F2»L3 ,RAP»L1 »LA»f P 1 *CFS*Sl *S3 
REAL Ll»L3.LA 
DO ?0 I« 7,9 

C EULER ANGLES IN RADIANS 
Yin » Yt 1 ) /RAD 
Ytl+Q) * Y(!+o>/RAD 
YII+1RI* Y( I + 1 8 1 /RAP 
C ANGULAR VELOCITIES IN RADIANS PEP SEC 
L « 1-3 

YCL I * YILI/RAD 
YIL+9) ■ YIL+9I/RAD 
70 YIL+1B1* YIL + 1B1 7RA>\ 

RFTIJRN 

FND 


Figure A36. Subroutine TOR AD Source Listing 



SU°ROUT I NE TOOEG 

C SUBROUT i,'4F TOOEG CONVERTS ANGLES AND ANGULAR VELOCITIES TO DEGREES 
COMMON/AAR/VI33I 

COMMON/AAH/C l *C3*F? »L3*RAD«L1 »LA*CE 1 *CF3» Si »S3 
REAL LltL3tL4 
C EULER ANGLES IN DEGREES 
DO ?00 I « T»o 
Y( n * Yt I >*RAO 
YII+R1 = Yl I ♦ RI*RA r ' 

VII+1AI* Y( I ♦lfl»*RAn 
C angular VELOCIT t e< in nrr-PFFS PER SEC 
L = 1-3 

YIL) * YCL I *RAD 
Y(L-*-OJ « Y(L+°I*RAn 
200 YIL-*-lf«l= Y«L+1B»*RAP 
RETURN 
END 


Figure A37. Subroutine TODEG Source Listing 


SUBROUTINE RANDUC YFL*SFED) 

SEED=AM0Dfl3l075»*SFED*3A3*>9738368.) 

YFL=SFED*.29l038304367F-ir) 

RFTURN 

END 


Figure A38. Subroutine HANDl) Source Listing 
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^tnppMTTNP OFOIVF C POL T * *&r #PHO *riUST ) 

rOMVO**/$A»/YI73l 

roMvm'/ur/^ ( > 

rnMvn^ /Aaj/**n^F 

r%| vrue; j qm OY(*»nt t 

->yvr^cTOM YOL'M*^) *oOLAt7r> 

ICAvrs^orr 

MP*^n 

r'O 70 T = 1 ,MC 
vCL^f T ) sY ( T ) 

^OL^f T I * rs f T ) 

TFrYm* r 0.r.f>)co TO 1A 
OYC n - # 01 *Y (II 

00 TO 15 

io ~vm * #oi 
IK ^nHTT v U c 
7n ^omttw’F 
wonc - 4 

CfOM = 1 # 

AO 70 

7* Cfftv = -c^rv 

opjt - tfr,« * n Y(,n 

Vfj) = vfjt + dcot 

r ^LL Ajpro<; 

ta^l r^UT* tru^Ti 

r«tl ro^TMPn 

r ali. foPfcc 

rA LL MC^c-MTC (PMA) 

OfFFOM 

Yfj) * Y (J) - proj 
TFtMON#F0#\ # > r O TO h* 
no -a* fs1 v NF 

01 ( T ) = ^(T ) 

-»*> 0(f> = fM7LA(T? 

OA TO 75 

4* 00 50 T = 1 »MP 

a i ff©* ncn - oicn 
1F(OTFP*#* P #0#0) n 0 TO 5* 

- c f 1 1 J) * % n 

on to 55 

*^o c f T # j ) * ATFFO/(7,*^Y( l) ) 
e«; rowTTWF 
5 0 rOMTTA»U p 

OO AS f * IfNF 
55 Of J) s on? Of f > 

70 rOMTfMMF 

rforwf - , «F,POl. T5 ) 

no 7M * t * NF 
o(n*ooLocn 
75 Yrn * YOLO(T) 

MOOFs f 5AV C 

pFTtrOM 

FNO 


Figure A39. Subroutine DERIVE Source Lasting 



SURROtlTINF F T OFN ( A • N »CHRONOS 1 
OIMFNSION AOO»lO) *RI90) »STORFI 2*0) 
1VAR(in» 

U.OV«*0*LI M !T 

L|WTT«L! v 'TT4-1 

TI*F(LIMIT)«CHR0N05 

cull mf5sfsmn»*»ni 

r«UL ORCAlt <*»A»R»".W1 


,TIMF(5I 




WlM>W 

TFILTWTT.FO.M r - n To 1 e 
of* 10 K*1» M 

10 5T0P r t U.OW-*-Xl*Rf*1 

orttiRM 

15 «TOP=o 

Of* ?0 t*i»5 

1F(MW( I) .r.T. w T0°1 »'Tf»P=VMft1 
?0 fOHTIMUR 

WRt TF (6*100) 

WR|TFI6»105> (TI**P« 1 )»I = 1 »*>» 
HRf T c (6*1101 
p»» 0(i)*10M(/‘»X» e ‘10.4 
rv^o( *0)*1 1"* »=X»fio.6) 

no 11 I=?t« 

11 IVARI n a 10W*1X*F10*4 
OO =0 IC=1»*T0P 


II*?»* 

I=TT-1 
JJ*60+II 
J=JJ-1 
LL»1’0+IT 
t -LL-1 
«M*1«n+! t 

MlsMM-1 

|Ftr,L F .*'»»<1 1 1PO TO 15 
«;T0PPtT1«^T0RFU n*10H 
IV ARC 1 >*10Ht /»1X»A10 
TVAR I ? ) *10H» IX* A1 0 
15 tF(r,LP« M V(?n^o to an 
CTORF ( J) sITOPFI JJ1 =10H 
!V«R(i) = tVAP.<a) = 10H*lx»iin 
40 IPIK.LE.MMC^nGO TO 45 
5TORP(I-1«STORFtUL1«inH 
IVAR(5)=IVAR I 5>=10H»1X*A10 
45 IF|*»IF*mm| 4> iGO TO 50 
eyOPF (»»i i=5T0RFC m R1*10m 
T'/»P(T 1sIV/*R(«1*10H»’X»»ln 

50 I F ( r . (_F « **** (5)1 Of* TO 55 
Rm=RUT)*lOH 
TVARIR>=10H*1X»A10 

iv«pnoi=ioH ,ix»#ini 

=0 PRINT lVAR*STORF(I 1 »STORF II 1 1 »STORF( J) »STOR c ( JJ1 »STORF M- 1 * 
1 5T0RFCLL I *STORF (N1 ) * STORE I NN) »R( I ) *RI 1 1 1 


Figure A40. Subroutine EIGEN Source Listing 



RFTURN 

100 F0RM«T(lHl*50x*11HP!r,FM\/ALnP«) 

10*> FORMAT* //**X*M*MTtMF =,Ff>,2*nX)) 

110 F0R^ATC/*^X*S(AHRFAL*9y* OH TMAG!NARY*&X) 1 
FNO 


Figure A40. Subroutine EIGEN Source Listing 
(Concluded) 



StJ9R0»*TINF HF*SrN(NtA*0) 

OIMrmMON hi n 
TAjTrr^P P*P*S ox f ^ 

TTf m # | pc TMPM 

*0*0+1 

**N*cw-l)*!n+t 

rx«w-Tn-io*i 

PX*M 

PM*1 

on 7i r * ?*K*tTn 

*'<*PX 

DMxOM^n 

oxaov+n 

JP*PM 

"*n # 

KsJP 

JX*J 

->0 T*A«SfMjn 

?*(T. !.*•"! '■’O TO ^S 

JCmjO 

JK = J 

n aT 

** TF(J*^F#MX> 00 TO **7 

,f0,fO + O 
OO TO *»0 

*** fC f ,J** F 0 # ri ao to a/, 

J*jr 

OO pspM f O)f 

Tan (O) 

MP1«A( J* 

MJ) = T 
*• J*J+1 

Osjr 
ao 

T*M.M 

M 

*<D)tT 

*>0 OsOf'' 

A4 TPf Af r) # ro # o # ) r.o TO 70 
jr«OM+n 
JK*v+1 

'** °zh<jr) 

T F( r .PO # 0.) oo TO 6* 

n a o*T 

»<rM*v+n 

j*n * jx+n 

on SO jM=JlfO f MM # 0 

AJM*M JM}- 0 ** C^ v ) 

rri nnc( ^jm> # L r .( Mjv) ) ) ) 

* f J*M =A J M 
to tCMsW+O 


HFS510P1 

HFSS100* 

MFSS1AP3 

MFSS1004 

HFSS100S 

MFSS1006 

MFSS1007 

HFSSlftOS 

HFSM"09 

MFSSlOll 

HFS51P1? 

HFSSini) 

HFSS1014 

MFSS101S 

HFSSJOU 

MPSSI 017 
HFSS1M8 
NFSS1019 
HFS5107P 
MF S S 1 0 ? 1 
MFS5102? 
HFSStO?3 

H*ssin?t 

HFSS1075 
HFSM 0?6 
NFSS1077 
HFSS107P 
HF SSI 079 
HFSS1030 
HFSMM1 
HFS5103? 
HFSS1033 
HFSS1P34 
HFSS1M5 
HFSS1036 
HFSSIM? 
HFS51098 
HFSS1A30 
Hfr<iM04n 
HFSS1041 
HFSS104? 
HFSS1A43 
HFSS1044 
MFSS1045 
HFSS1A46 
HPSS1047 
HFSS1048 
HFSS 1049 
HFSS1080 
MFSS10S1 

Ajvo0 t MFSS10S? 

HFSS1033 

HFSS1034 


Figure A41. Subroutine HESSEN Source Listing 
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jor 

HFSM065 


'Vl 60 P«PM # PX 

HFSS! 056 


AP*A(P)+ft«At J) 

MFSMP57 


IPeAO<;(ADJ.LF.(.lP-<»*An«;(A(P) ) ) ) ApaO, 

MPSF] 0*A 


A ( P ) aAP 

HFSM05O 

60 

J»J*1 

HFFS1060 

6* 

JKajr-M 

HF«5Fin61 


jr x jr+rv 

MFSF106? 


TF{ JK.l_F.MtO '■-O TO 4* 

MPV<1«6* 

70 

PONT ?•>>!)* 

He , 'K1064 

71 

COWTTMII* 



OFTtIPM 

MCSF1065 


FND 

HFSM066 


Figure A41. Subroutine HESSEN Source Listing 
(Concluded) 
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c»fooOHTl**F OPCALUr* * 

nturvcfn*' A ( 1 ) i° ( 1 ) 

T My corn «> 

M « M T M 
AVM = 1 # 

ATT s . 7 P—7 
Arc * # ir-to 
tyc*> * o 
'• a n 
MO - n+i 
MM - M + ( M_ 1 ) *C 

VM s 

Vt a 

TF(»'^0 # 1 ) r . 0 TO if 

T r C M »t r # , > orTI ,0K * 

TF(»* # rr # -»> or yr *>* 

k nc | t i = act>aoc ( * (mm ) ) 

^I.T 1 ' = ^MAXM np LTA f Af*“> 

\rr a c A (WV +1 1 ) 

TC ( *or # rr).r. ) or Tf> if 
tcprr/.^^TiTf) r.O TO 7* 

T ( tT^D.OT.^ \ on TO i f 
TF( rn y^ '> c 

T «. M . s* + ? 

V r 

!F( ap«(X WL c . # 1 r -*> X*% 

o f > =x 

^ f *M - ' « 

* T 1/ * »iT*UM4,1 

t r f 5 { M UP^,A # ) r.o to ic 

V a - P("-i ) /y 

on re *>i 

ic ^rvTf'i'c 
-»*» T ▼ TO - * 

a^t = # ir-7 

*■ = M.1 

MM — MM 

m»* a fc, L 

*k‘| s 

! f l M #r T i^' or t n ^ * 

TC(M # rn # -)l rp TO ->R 

tC(^ # co#i ) on yn if 

‘'sV/1 
»~Tt !*>•• 

*>r r - .*#( A < W»M*A ( H(t* ii 

c>m - ARC ( f ( mm ) ~!\ ( mw ) i 
MM * A»<(A(NW)HAP«(MN^n 

f c ( n/» m # lf* a r t#*am ) OAA»a% 

n.^A’an 

^ r A f mm+i i * A ( VM.i ) 

t s n *v + r 

Tree A"MT) .LF. Apt#*vm) # 0P* ( a°MT 1 •L c ’»A'“T# n *'M ) T*% 

O a COOT ( Aoc; ( T ) ) 


ORCA 1001 
0RCA100? 
ORCAt 003 
OPC A 1004 
0PCAin*5 
opcaioo* 

0PCA1AA7 

OPCAIOOR 

OPCAIOOQ 

ORCA101* 

ORCAl^n 

OQCA101? 

ORCAini? 

0PCA1014 

ORCA'.ri* 

OPCAlOlf 

oocaioi7 
ORCAI^l R 
0RCA1M° 
0PCA1O70 
r»C A!0?1 
ORCAIO?? 
0PCA107* 
0RCA10?4 
ORCAin?* 
OPCAlo?ft 
noCAin?7 
OPCpinro 
norAir?p 
0PCA103O 
0PCA10P1 

npCAl o3f 

0RCAir3A 

0RCAI039 

0PCA104P 

0PCA1O41 

nPCAlOA? 
0PCA!04^ 
OPCA1044 
0PCA1045 
0PCA1046 
OPCA1047 
0PCA1 04 A 
ORCA 1 040 
ORCA 1 Of r 
0RCA1051 
ORCA1052 
0PCA10A3 
0PCA1 0*4 
OPCA 1 off 
OPCAinAf 
OPC A 1 0 * A 


Figure A42. Subroutine QRCALL Source Listing 


169 



no TO *0 

TF(T*rF,n*| 0-0 TO *n 

*' * »U’ 

o(M-l ) = n 

O t M ) * r 
7 7 M « N-! 

MM S MM 
*‘M 3 M|_ 

M L = VL- W * 
no ta 17 

-»n »' a 

D(M.i ) s o + r 
O(M) s n # 

* * MfM-M+l 

** a '♦+? 

OfV-1) = n_r 
O(M) = f) m 

no to *>7 

*o rr ( t. '■•=•. n. i r-n to a* 

O ( M+c J = Q 

-CM.*) = O 

^{‘U7j z r 
ofM+B] - _r 

no to 70 
ah x a » 4 .r 

V * n_r 

= n. 

'JfM+O) - 0, 

9(M45) = X 
9(M+7) a v 

|F( AOnfXl.rT.APMY) 1 00 TO 7r 

r?fM+R) . y 
= X 

7 ^ TF( TTt-R.t.r^n ) r.^ 

y , AQC{9(M+* )-9 (M4.1 ) ) + A r >< (P (f/4-ft)- 1 ' f ” + ? ) ) 

Arc * a«mp< m *m> +AOMM M +6n 

TF( Arr*r,T # i # > x*x/ACr 

Y a AOS(P<M+7 t-P<M+M )+A«MR(M+B)-P<? 4 *M 1 

arr = i°«;(p(^4.7 ) i ♦^mocumi 

TF< Arr.r,T # t. ) y*y /A cr 
Arc * a*m mnl*i i i 

OCLTA * AMAX1 (OFLTAt ( ArT*AP.0( A(NM)1 H 
TFf Arr t r,T«OPLTA ) ro TO *n 
TF(TT-R*0T # 4 ) r.o TO 

!F( ( X . L p ♦ ATT ) #AA»o,(Y#L r #AnT) ) no TO 7* 

on ]P\ jTrp.rT.^nn) ( To r>n r 

TFffX.rT#.* uav%(y # at m « 1) no to no 

X * m +5 

TFfY.rT..* ) no to j?o 

TF ( X.CT ♦ . A ) 00 TO no 

puo « R(v+M«P('47^q(^M^( , >®) 

CfOM a 0(*>+« j+P(M+?) 

no roMTiwtiF 

ANN ■ A ( NN ) 


0PCA10** 

OPCA1A60 

OPCA1061 

0PCA106? 

OPCAlOf^ 

0PCA1064 

0PCA106S 

opcaio** 

0 »CA!o 67 

O&rAlOAA 

ODFAlOftO 

OPC At 070 
OOCA1A71 
OPCA) 077 

0PCA7074 
0PCA1O76 
0PCA1077 
OPrA107« 
OPCA\07P 
opr a 1 n »o 
nor A 1 Ofll 
oorftl 007 

OPfATOfto 

OprAinfU 

oprAio**? 

noCAi oftf 
0PCA1O47 
Opr A 1 n a 0 
o CA70P9 
jPCAIoor 
OPCA J OP 1 

OPfAl f\Q? 

OPCA 1 00 0 
OPOAl A04 
OPCA 1 00^ 
OPCA 1 00ft 
0PCA10O7 
OPCA1000 
OPCA 1 OOP 
OPCA 11 Op 

OPCA 1101 
0PCA110? 
oPCAiim 
0 PCA 1 1 04 
0PCA1 1 A4 
OPCAIIO^ 
OPC A1 1 07 
OPCAUO* 
OPCAinP 

oPCAino 

opcAnn 

OPCAtn? 

0PCAU13 


Figure A42, Subroutine QRGALL Source Lieting 
(Continued) 
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Tfutpr.lf.isi r, o to in? 

!F(PMO,NP.n* > no to in? 

lF(*fr,wa # *|s,n # ) no TO 109 

?mo ■ .nom 

ato*'a ■ .n? 

10? CALL OR(N.A»PHO.«iIG u A.n.OP|.TAl 
*» > A»r>l AtNN) 1 
ANN ■ AP5( ANN-AtNNl 1 
IF( °.PT » ACT) ANN ■ ANN/° 

TTPO « ITPR+1 

IF < IT e R»FO,?5 1 a<-t«.!P-« 

OO T»l*4 
K » *♦! 

10* R(V) ■ RTK+A) 
po TO 1 * 

1 1 n r • »•+■> 

1?n RMO * R{P )*R ( p > 

MO**a « Rt < ) (p i 
PO TO 100 
1*0 RHO «• 0, 

«IPMA * 0. 
po TO 100 
*nn covTiN'in 

MRt TF lf>*?10) M1W 

?10 PORMATJ/AH all .U»?4H FlfiFNWALHFS 
'4«0 

orTiio* 1 

p«n 


0RCA111A 
OOpAI 115 

ORCA1 l 1ft 
09PA1 117 
OorAi 1 1 A 

0RCA\11« 
0RCA1 1 ?o 
OPCA1 1 ?1 
09C A 1 1 ?? 
ORPA1 1 ?5 
ORCA!l?A 
ORPA1 1 ?5 
0»CA11?A 
0RCA1 1?7 
ORPA 1 1 ?fl 
OPPA1 1 ?« 
OOCA1 I7f 
ORP A1 1 *1 
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Figure A43. Subroutine QR Source Listing 
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